Exploiting the Einstein Telescope to solve the Hubble tension

The detection of Gravitational Waves (GWs) from the coalescence of merging Binary Black Holes (BBH) and Binary Neutron Stars (BNS) [1, 2] opened a new window to test the General Relativity, relativistic astrophysics, and cosmology [2, 3, 4]. As it is well known, the GWs bring direct information on the luminosity distance of sources and, therefore, they can be used as rulers to measure distances in the Universe. Indeed, they are usually called standard sirens [5, 6], and are fully complementary to the standard candles, such as Cepheids and Supernovae Type Ia (SNeIa) among the others, that are instead based on the detection of their electromagnetic emission and need to be calibrated on closer sources in order to get a measure of their luminosity distance. Although GWs offer an alternative method to obtain distances in cosmology and are not affected by calibration problems, they are not free of issues. Indeed, the GWs waveform encodes both information on the systems, such as masses, spin, and inclination angle among others, and information on the cosmology such as the distance. Furthermore, they encode information on a given theory of gravity and, potentially, can probe it [7, 8, 9]. However, information on masses and redshift is completely degenerate, and the only way to break such degeneracy is to have prior information on the redshift from an electromagnetic counterpart. There are several ways to get accurate information on the redshift. For instance, one can assign to the GWs source the redshift of the host galaxy [5, 6, 10] or, alternatively, looking at the electromagnetic emission following the GWs such as short Gamma-Ray Burst (GRB) [11, 2] and kilonovae [2]. However, the host galaxy can be accurately detected only up to redshift below one [12], and kilonovae will be detected up to redshift $z\sim 1$ [13]. On the contrary, short GRBs may be detected using forthcoming satellites, such as Transient High Energy Sources and Early Universe Surveyor (THESEUS), up to redshift $z\sim 8$ [14, 15, 16]. Such high redshift detections can allow also to deeply test the cosmological evolution [17, 18, 19, 20]. Furthermore, a complementary avenue to obtain the redshift information is represented by the observation of tidal deformation in BNS mergers [21, 22]. Indeed, it may supply redshift estimation with an accuracy ranging from $8$% to $40$ % depending on the choice of the equation of state (EoS).

Nowadays, the LIGO/Virgo/KAGRA collaborations have explored several ways to constrain cosmological models from GWs events. A turning point was the event GW170817, i.e. the first merger of a BNS with the simultaneous detection of the GRB 170817A [2], yielding to the first estimation of the Hubble constant with GWs, $H_{0}=70_{-8}^{+12}$ km s^-1 Mpc^-1 at 68% of confidence level [23]. Afterward, the LIGO/Virgo/KAGRA collaborations have explored the possibility to constrain $H_{0}$ by analyzing the population distribution of BBH mergers, and searching for the host galaxy identification [24]. However, these new measurements of the Hubble constant are in agreement with both the late-time and the early-time measurements of $H_{0}$ and, therefore, do not help to solve the so-called Hubble tension [25, 26, 27, 28, 29] which may, however, be alleviated to a 2.1$\sigma$ tension if the errors on the Hubble constant were underestimated [30]. Other examples of solutions to the Hubble tension may rely on modifications of fundamental laws [31, 32, 33]. Nevertheless, the next generation of GW detectors, e.g. the Einstein Telescope (ET), can strongly improve the accuracy on the Hubble constant reducing it below 1% [34, 35], and promise to offer a solution to such a tension pointing out its correct value. Therefore, there is an important need to study also the theoretical framework related to the Hubble tension. Let us remember that the Hubble tension is $4.2\sigma$ discrepancy between the measurements of $H_{0}$ obtained by fitting the CMB power spectra [36], and using the standard candles such as Cepheids [37] both in the framework of the concordance cosmological model, also known as $\Lambda$ Cold Dark Matter ($\Lambda$CDM) model.

Since the nature of Dark Energy (DE) is still a puzzle (for a comprehensive review see Ref. [26]), there are many attempts to explain the $H_{0}$ tension modifying the DE EoS [38, 39, 25, 26, 40, 41], or the underlying theory of gravity [42, 43, 26, 44]. Several analyses have been carried out to measure the Hubble constant and DE EoS with GWs making use of different techniques such as the identification of the electromagnetic counterpart of GW sources [45, 46, 47], the cross-correlation between GW sources and galaxies [48], the statistical host identification techniques with galaxies [12, 49, 50, 51, 52], the hierarchical inference without galaxy surveys [53, 54, 55], the study of lensed events [56, 57] and the analysis of NS EoS [58, 59, 60].

It is worth noting that the $\Lambda$CDM model has its own issues such as, for instance, the well-known small-scale problem of the CDM paradigm and the extremely small value of the cosmological constant compared to the expectation from quantum field theory [61]. Therefore, alternative approaches have been considered to overcome both the lack of experimental detection of the dark ingredients and the small-scale issues of CDM as well. For instance, the modification of the underlying theory of gravity can explain the acceleration of the Universe by means of extra degrees of freedom arising from higher-order curvature invariants or extra scalar fields [62, 63]. Alternatively, models of non-homogeneous universes or violations of the Copernican principle can also account for the accelerated expansion [64, 65]. Nevertheless, the debate on whether there is or not a real need to shift from the $\Lambda$CDM to a more complicated model is still under debate [66].

Here, we will focus on a set of models which modify the DE EoS in view of solving the Hubble tension [25, 26]. For each model, we will predict the accuracy down to which ET will be able to detect departures from the $\Lambda$CDM offering in such a way a theoretical framework of DE capable of solving the Hubble tension. In Sec. II, we will briefly introduce the DE models. In Sec. III, we will summarize the procedure adopted to build mock data that will mimic the ET observations of the luminosity distance. In Sec. IV, we will give details of our statistical analysis, while in Sec. V we will show results. Finally, in Sec. VI we will give our final discussion and conclusions.

We focus on four DE models which may help in solving the Hubble tension [25, 26, 67, 68], and differ from each other in the way they affect the DE EoS leading to a modification of the luminosity distance. In the context of General Relativity, and of the Friedman–Lemaître–Robertson–Walker (FLRW) cosmology, the luminosity distance is defined as

$$d_{L}(z)=\begin{cases}\frac{c(1+z)}{H_{0}}\frac{1}{\sqrt{\Omega_{k}}}\sinh{\left[\sqrt{\Omega_{k}}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\right]}\quad&\mbox{for}\ \Omega_{k}>0,\\ \frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\qquad&\mbox{for}\ \Omega_{k}=0,\\ \frac{c(1+z)}{H_{0}}\frac{1}{\sqrt{|\Omega_{k}|}}\sin{\left[\sqrt{|\Omega_{k}|}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\right]}\quad&\mbox{for}\ \Omega_{k}<0,\end{cases}$$

(1)

where $z$ is the redshift, $c$ is the speed of light, $H_{0}$ is the Hubble constant, $\Omega_{k},0$ is the value of the curvature parameter at $z=0$, and

$$E^{2}(z)=\Omega_{m,0}(1+z)^{3}+\Omega_{k,0}(1+z)^{2}+\Omega_{DE}(z)\,,$$

(2)

where $\Omega_{m,0}$ is the value of the matter density parameter at $z=0$, and $\Omega_{DE}(z)$ is the DE density parameter as a function of the redshift. In the $\Lambda$CDM model, $\Omega_{DE}(z)=\Omega_{\Lambda,0}$. Since Eq. (2) is strictly related to the Friedman equations and to the DE EoS, changing the DE model leads to different expressions of the function $E(z)$ [69]. To study the proprieties of DE component in a general framework, one can consider the ratio of the DE pressure to its energy density as a function of redshift [70, 71]:

$$\omega_{DE}(z)=\frac{p_{DE}(z)}{\rho_{DE}(z)}\,.$$

(3)

Again, $\Lambda$CDM model is recovered for $\omega_{DE}(z)=-1$. In the next subsections, we will discuss four DE models whose modifications of the Eq. (2) may serve to solve the Hubble tension. Moreover, we will also discuss the limits in which such models recover the $\Lambda$CDM cosmology.

Here we briefly summarize the procedure adopted to build up the mock catalogs. We closely follow the recipe given in [91], and assign the redshift to the GWs sources extracting it from the following redshift probability distribution [92, 93, 39]

$$p(z)=\mathcal{N}\frac{R_{m}(z)}{1+z}\frac{dV(z)}{dz},$$

(19)

where $\mathcal{N}$ is a normalization factor, $dV(z)/dz$ is the comoving volume element, and $R_{m}(z)$ is the merger rate per unit of volume in the source frame. The latter takes the form [94, 95, 96]

$$R_{m}(z)=R_{m,0}\int_{t_{min}}^{t_{max}}R_{f}[t(z)-t_{d}]P(t_{d})dt_{d}\ ,$$

(20)

where $R_{f}[t(z)-t_{d}]$ is the Star Formation Rate (SFR), and $P(t_{d})$ is the time delay distribution. We assume, for the SFR, the model proposed in [97] and for the time delay distribution a power law functional form, $P(t_{d})\propto t_{d}^{-1}$, as suggested by population synthesis models [98, 99, 100, 101, 102]. Nevertheless, it is worth noticing that setting the SFR and the time delay distribution to other models do not affect the accuracy of final results (for more details we refer to Sec. 5.3 in [91]). We integrate Eq. (20) between a minimum time delay of $20$ Myr and the maximum fixed to the Hubble time. Furthermore, the quantity $R_{m,0}$ is the normalization of the merger rate at $z=0$. We, therefore, set it to the best-fit value obtained by the LIGO/Virgo/KAGRA collaboration: $R_{m}(z=0)=105.5^{+190.2}_{-83.9}\ \mbox{Gpc}^{-3}\mbox{yr}^{-1}$ [103]. Once the redshift is extracted from the probability distribution in Eq. (19), we can assign a fiducial luminosity distance, $d_{L}^{fid}(z)$, based on our fiducial cosmological model in Eq. (18).

The ET will have three independent interferometers and, hence, the combined SNR is $\rho=\left(\sum\limits_{i=1}^{3}(\rho_{(i)})^{2}\right)^{1/2}$. The SNR of the single interferometer, $\rho_{(i)}$, in the ideal case of Gaussian noise, is:

$$\rho^{2}_{(i)}=4\int_{f_{\rm lower}}^{f_{\rm upper}}\frac{|F_{+,i}\tilde{h}(f)_{+}+F_{\times,i}\tilde{h}(f)_{\times}|}{S_{h,i}(f)}df.$$

(21)

In the previous definition, $S_{h,i}(f)$ is the one-sided noise power spectral density of $i$-th interferometer, $\tilde{h}_{+}$ and $\tilde{h}_{\times}$ are the GW strain amplitudes of $+$ and $\times$ polarizations, and $F_{+,i}(\psi,\theta,\phi)$ and $F_{\times,i}(\psi,\theta,\phi)$ are the so-called beam pattern functions [104]. The whole sensitivity function¹¹1The latest power spectral density $S_{h}(f)$ can be downloaded at https://apps.et-gw.eu/tds/?content=3&r=14065. $S_{h}(f)$ is depicted in in Fig. 2 . To integrate Eq. (21), we set a lower cutoff, $f_{lower}$, at $1$ Hz [105] and the upper one to $f_{upper}=\frac{c^{3}}{(6\sqrt{6}\pi M_{obs})G}$.

We can compute the total number of observable BNS mergers, $N$, from the equation

$$N=T_{obs}\ \Theta\int_{0}^{10}\frac{R_{m}(z)}{1+z}\frac{dV(z)}{dz}dz\,,$$

(22)

where $\Theta$ is the duty cycle and $T_{obs}$ is the observation time. In order to generate the catalog and select the events above a certain threshold of the signal-to-noise ratio (SNR) $\geq 9$, we assume an isotropic distribution for the sky angles $\theta$ and $\phi$, and an uniform distribution for the orientation angle $\cos i$ and the polarization $\psi$. Moreover, to generate the synthetic signal self-consistently with our choice of $R_{m,0}$, we follow the LIGO/Virgo/KAGRA collaboration and set a uniform NS mass range in the interval $[1,2.5]\ M_{\odot}$. Thus, we obtain a rate of $\sim 3\times 10^{4}$ events per year, assuming a duty cycle of 80%. In our analysis, we will consider 1, 5, and 10 years of observations.

Once the fiducial luminosity distances are generated, we add a Gaussian noise component, $\mathcal{N}(d_{L}^{fid},\sigma_{d_{L}})$, to them in order generate our mock observations. The variance $\sigma^{2}_{d_{L}}$ includes the contributions due to the instrument, $\sigma^{2}_{inst}$, the lensing, $\sigma^{2}_{lens}$, and the peculiar velocity of the host galaxy, $\sigma^{2}_{pec}$. Therefore, the total variance will be:

$$\sigma_{d_{L}}^{2}={\sigma_{inst}^{2}+\sigma_{lens}^{2}+\sigma_{pec}^{2}}\,.$$

(23)

The contribution due to the instrumental noise component $\sigma_{inst}$ is [106, 107]

$$\sigma_{inst}=\frac{2}{\rho}d_{L}(z),$$

(24)

where factor two accounts for the degeneration between $\rho$ and the inclination angle, which may differ for each event. The contribution due to the weak lensing distortions, $\sigma_{lens}$, is given by [108, 109]

$$\sigma_{lens}=0.066\left(\frac{1-(1+z)^{-0.25}}{0.25}\right)^{1.8}d_{L}(z)F_{delens}(z),$$

(25)

where $F_{delens}(z)=1-\frac{0.3}{\pi/2}\arctan{\frac{z}{z_{*}}}$, with $z_{*}=0.073$ [109]. The latter factor takes into account the possibility to reduce the uncertainty due to weak lensing with the future detectors such as the Extremely Large Telescope [110]. Finally, $\sigma_{pec}$ is related to the peculiar velocities [111, 112, 113] and can be approximated with the following fitting formula [114]

$$\sigma_{pec}=\left[1+\frac{c(1+z)^{2}}{H(z)d_{L}(z)}\right]\frac{\sqrt{\langle v^{2}\rangle}}{c}d_{L}(z)\,,$$

(26)

where we set the averaged peculiar velocity $\langle v^{2}\rangle$ to $500$ km/s, in agreement with the observed values in galaxy catalogs [115]. To build such a mock data set we are relying on the somehow simplistic assumption that there are no correlations to take into account with the source parameters of each system. Of course, in reality, this will be not the case. However, following [106, 107, 93, 39, 38, 67, 116, 117, 118], we are assuming that the luminosity distance is measured with a statistical error that depends on the observational features of the Einstein Telescope. Moreover, it is worth remarking that our estimation of the error bars is, on average, a factor $1.7$ larger than the one expected for the Einstein Telescope [34] providing, therefore, a conservative estimation of the luminosity distance.

We carry out a Monte Carlo Markov Chain (MCMC) analysis to estimate the accuracy down to which each DE model parameter (that depends on the specific choice of the DE EoS) can be constrained with the future observation from ET. Our mock data are built using the flat-$\Lambda$CDM cosmology in Eqs. (2) and (18) as our fiducial model. Then, we expect the posterior distributions of the parameters of the DE models introduced in Sec. II to be centered around the fiducial model. Therefore, the error on the model parameters will indicate the accuracy that we will be able to reach with ET. To this aim, we will run our MCMC pipeline on both bright and dark sirens, i.e. events whose electromagnetic counterpart has been and has not been detected, respectively, and we will point out the main differences in the results.

Our MCMC is based on the emcee package [131], and will employ the likelihood of all GWs events defined as the product of the single event likelihood, $p(\textbf{d}|\bm{\lambda})=\prod_{i=1}^{N}p(d_{i}|\bm{\lambda})$. Here, $\bm{\lambda}$ are the cosmological parameters of interest for the specific model, and d$\equiv\{d_{i}\}_{i=1}^{N}$ is the mock dataset with $N$ equal to the number of observations. In order to write down the single event likelihood, one must distinguish between the run with the bright and dark sirens [91]. When using bright sirens, the redshift information is assumed known from the detection of an electromagnetic counterpart which is, in our case, a short GRB. In such a case, the single event likelihood can be written as [132, 54]

$$p(d_{i}|\bm{\lambda})=\frac{\int p(d_{i}|D_{L})p_{pop}(D_{L}|z_{i},\bm{\lambda})dD_{L}}{\int p_{det}(D_{L})p_{pop}(D_{L}|z_{i},\bm{\lambda})dD_{L}}\,,$$

(27)

where $p_{pop}(D_{L}|\bm{\lambda})=\delta(D_{L}-d_{L}^{th}(z_{i},\bm{\lambda}))$ [133]. In the Eq. (27), the denominator is a normalization factor that takes into account the selection effects [132, 134].

Instead, when using dark sirens, we assume that the redshift distribution of the BNS population is known, and marginalized over the redshift [55, 54]. In such a case, the probability of detecting an event $d_{i}$ in a specified cosmological model is given by

	$\displaystyle p(d_{i}|\bm{\lambda})$	$\displaystyle=\int_{0}^{z_{max}}p(d_{i},z_{i}|\bm{\lambda})dz_{i}$		(28)
		$\displaystyle=\int_{0}^{z_{max}}p(d_{i}|d_{L}^{th}(z_{i},\bm{\lambda}))p_{obs}(z_{i}|\bm{\lambda})dz_{i}\,,$

where the probability prior distribution of the redshift, $p_{obs}(z_{i}|\bm{\lambda})$, is obtained from the observed events and already includes detector selection effects [55].

Finally, in our analysis, we neglect (i) the contribution of the spin of the source to the amplitude of the signal [135, 136], (ii) assume a flat uniform prior on the cosmological parameters of interest as reported in Table 1, and (iii) we neglect the impact of the merger rate and of the time-delay distribution for being negligible within this dataset (we refer the reader to Sec. 5.3 in [91]).

We carried out a MCMC run for each mock catalog and for each DE model. We built six mock catalogs: three of them reporting all the GW events at one, five, and ten years of observations but without having prior redshift information (dark sirens); and the other three containing those GW events with a detected electromagnetic counterpart and, therefore, having prior redshift information (bright sirens) from a X-ray telescope, such as THESEUS. We use these mock catalogs to constrain the four DE models mentioned in Sec. II. All results of our run are reported in Table 2. For each DE model we also show the corner plot of the posterior distributions of the cosmological parameters of interest, see Figs. 4, 5, 6, 7, and 8. In each Figure, we report on the left and right sides the posterior distribution obtained from the bright and dark sirens, respectively. In each contour plot, we depict with green, orange, and blue histograms and filled areas the results from one, five, and ten years of observation, respectively. The different level of the transparency of the contours corresponding to a specific total number of years of observations depicts the 68%, 95%, and 99% CL from the darkest to the lightest color, respectively. Finally, the vertical dashed red line indicates the value of the fiducial cosmological parameters. Let us now discuss in detail the results for each DE model comparing them with the current observational results. In the discussion, we will always refer to the results obtained after ten years of observations.

$\omega$CDM: we report the results in Fig. 4. As reported in Table 2, we may constrain the cosmological parameters with an accuracy of $[\sigma_{H_{0}},\sigma_{\Omega_{k,0}},\sigma_{\Omega_{\Lambda,0}},\sigma_{\omega_{DE}}]=[1.02,0.18,0.18,0.93]$ and $[0.06,0.02,0.03,0.10]$ in the case of bright and dark sirens, respectively. These results correspond to a relative accuracy of the Hubble constant and the $\omega_{DE}$ of $[\sim 1.5\%,\sim 69\%]$ and $[\sim 0.1\%,\sim 10\%]$, in the case of bright and dark sirens, respectively. As a reference, the current accuracy on $H_{0}$ is at a level of 1%, and on $\omega_{DE}$ is at a level of 3% [73]. However, it is important to remark that those bounds were obtained by using not only luminosity distances but also the distance prior data obtained by Planck satellite [73]. Therefore, we show that by using bright sirens, ET will bound the Hubble constant with the same level of accuracy, but it will be capable of improving it by one order of magnitude using dark sirens. Nevertheless, ET alone will not be capable of improving the accuracy on the $\omega_{DE}$ parameter not even with the dark sirens catalog.

Interacting Dark Energy ($\omega_{DE}$-fixed): we report the posterior distributions in Fig. 5, and the best fit values of the cosmological parameters in Table 2. The 68% uncertainties are: $[\sigma_{H_{0}},\sigma_{\Omega_{m,0}},\sigma_{\xi}]=[0.95,0.13,0.88]$ and $[0.05,0.01,0.06]$ in the case of bright and dark sirens, respectively. These results translate in accuracy on the Hubble constant of $\sim 1.4\%$ and $\sim 0.1\%$, which is always better than current constraints shown in [79, 80, 81] by a factor of $\sim 2.4$ and $\sim 46$. On the contrary, the accuracy on the parameter $\xi$ is improved only with dark sirens by a factor of $\sim 3.3$.

Interacting Dark Energy ($\omega_{DE}$-variable): the results of the MCMC algorithm are shown in Fig. 6 and listed in Table 2. The cosmological parameters are constrained with an accuracy of $[\sigma_{H_{0}},\sigma_{\Omega_{m,0}},\sigma_{\omega_{DE}},\sigma_{\xi}]=[1.19,0.21,0.53,1.5]$ and $[0.04,0.16,0.08,0.17]$ which corresponds to a relative accuracy on the Hubble constant and the $\omega_{DE}$ of $[\sim 1.7\%,\sim 38\%]$ and $[\sim 0.1\%,\sim 9\%]$, in the case of bright and dark sirens, respectively. Using only bright sirens, the uncertainty and the relative accuracy on $H_{0}$, $\omega_{DE}$ and $\xi$ are not comparable with the ones obtained in [79, 80, 81]. Nevertheless, when we use dark sirens the accuracy on $H_{0}$ improves of a factor $\sim 27.5$ while the constraints on $\omega_{DE}$ and $\xi$ do not still improve results from [79, 80, 81]. It is worth noticing that previous analyses are based on multiple datasets, such as CMB and Cepheids, while we are only focusing on studying the capability of ET.

Emergent Dark Energy: the results are reported in Fig. 7 and in Table 2. We constrain the cosmological parameters with the following accuracy: $[\sigma_{H_{0}},\sigma_{\Omega_{m,0}},\sigma_{\Delta}]=[1.43,0.06,0.86]$ and $[0.03,0.002,0.01]$ in the case of bright and dark sirens, respectively. These results provide us with a relative accuracy on $H_{0}$ of $\sim 2.1\%$ and $\sim 0.04\%$. Using bright sirens allows us to obtain bounds on the Hubble constant comparable with current constraints shown in [86]. Instead, using dark sirens, we improve such a constraint by a factor $\sim 46$. The parameter $\Delta$ is constrained with a better accuracy only when dark sirens are taken into account, improving the bounds in [86] of a factor $\sim 40$.

Time-Varying Gravitational Constant: we report the posterior distributions in Fig. 8, while the constraints on the cosmological parameters are reported in Table 2. We show that, in the framework of a time-varying gravitational constant, ET will be capable of bounding the cosmological parameters with an accuracy of $[\sigma_{H_{0}},\sigma_{\Omega_{m,0}},\sigma_{\delta_{G}}]=[0.86,0.08,0.44]$ and $[0.04,0.01,0.02]$, in the case of bright and dark sirens, respectively. Hence, the predicted relative accuracy on the Hubble constant is $\sim 1.2\%$ and $\sim 0.06\%$. Using the CMB + BAO + SN +$H_{0}$ dataset, the bounds on the Hubble constant are currently at the level of $\sim 1.5\%$, while the accuracy on $\delta_{G}$ is of the order $0.002$ [73]. Therefore, while using dark sirens ET will be capable of improving more than one order of magnitude of the relative error on $H_{0}$, it will not be capable alone of improving the bounds on $\delta_{G}$.

The Hubble tension is one of the most important issues of modern cosmology [25, 26]. It is still not clear whether the solution to such a tension regards more the observational and statistical sector than the theoretical one with the possibility of some "new physics". Most of the solutions proposed up to date are focused on extending the $\Lambda$CDM model [39, 25, 26], and on changing the underlying theory of gravity [42, 43, 26, 44, 69]. Nowadays, this tension is established to be at $4.2\sigma$ and arises from a discrepancy in the value of the Hubble constant obtained from late-time observations, such as Cepheids, SNeIa, and BAO among the others, and the observation of the CMB power spectrum at early-time. A dataset complementary to the usual late-time observations is represented by the estimation of the luminosity distance from the GWs. Since the latter is not based on the measurements of the photon flux, they must not be calibrated on the closer electromagnetic sources, such as Cepheids and SNeIa. Therefore, they represent a potential way to solve the Hubble tension and may identify its cause whether it is related to the observations at the late or early time, or to the theoretical limitations of the $\Lambda$CDM model. To this aim, the LIGO/Virgo/KAGRA collaborations have constrained the Hubble constant with GWs to be $H_{0}=70_{-8}^{+12}$ km s^-1 Mpc^-1 at 68% CL [23]. However, the accuracy reached is not enough to provide a definitive answer.

LISA and 3G detectors such as ET will provide GW standard sirens at high redshifts, allowing us to test different dark energy models and, possibly, the dark energy equation of state. The population of BNS mergers detected out to redshift equal to 3 by 3G detectors, and the population of black hole binaries mergers that will be detected up to redshift 10 by LISA, will allow to probe the cosmological principle, to map the dark matter distribution, to probe the equation of state of dark energy [137, 138]. The 3G GW detector ET promises to constrain the Hubble constant with sub-percent accuracy [34], offering a possible solution to the Hubble tension. Therefore, we have forecast the accuracy down to which ET may bound the cosmological parameters of four DE models which have the potential to solve the Hubble tension [26]. Namely, we investigate the non-flat $\omega$CDM, the interacting dark energy, the emergent dark energy, and the time-varying gravitational constant models. We have predicted the luminosity distance expected in those models varying the cosmological parameters, and fit it to the mock data built to mimic the expected rate of observations and accuracy of ET, whose construction was explained in Sec. III. Our fitting procedure is based on the MCMC algorithm explained in Sec. IV. The results are reported in Table 2, and we also show the posterior distributions of the cosmological parameters of interest in Figs. 4, 5, 6, 7, and 8.

Our results clearly indicate that ET will be capable of reaching an accuracy of $\sim 1\%$ with bright sirens, and go below $\sim 0.1\%$ with dark sirens, independently by the theoretical framework used in the statistical analysis. This accuracy will be adequate to solve the Hubble tension. Nevertheless, ET alone will not always be capable of improving current constraints on the additional cosmological parameters that depend on the specific choice of DE model. For instance, in the non-flat $\omega$CDM and interacting DE models, the parameters $\omega_{DE}$ and $\xi$ will be constrained with an accuracy worse than current bounds [79, 80, 81, 73]. In the case of the time-varying gravitational constant model, the accuracy reached by ET will be still one order of magnitude higher than current constraints [73]. On the contrary, in the emergent DE model, we show that ET will be also able to improve the bounds on the additional cosmological parameter $\Delta$ by a factor of 40 with respect to current analysis [86]. These results show the huge capability of ET to solve the Hubble tension independently by the theoretical framework chosen but also point out that, to strongly constrain the DE models we have considered, ET will need to be complemented with other datasets.

For instance, LISA will be able to obtain an accuracy on the dark energy cosmological parameters similar to the one we forecast for ET [109]. On the other hand, LSST will provide weak lensing and BAO observations to set stringent constraints on the equation of state of dark energy [139]. Indeed, using weak lensing and BAO separately, LSST will achieve accuracy on the curvature parameter of $10^{-3}$ that can be strongly improved by a joint analysis of these two cosmological probes [140]. This represents an improvement in the accuracy of the curvature parameter $\Omega_{k,0}$ of three orders of magnitude with respect to our results. This is due to the nature of the cosmological observations provided by LSST that can probe curvature much better than the $d_{L}$-redshift relation. Finally, another interesting possibility will come from SKA which will allow testing cosmology up to redshift $z=6$ using neutral hydrogen intensity mapping [141]. Indeed, it will allow us to bound the dark energy equation of state with an accuracy of $\sim 0.34\%$ [142].

MC, DV, and SC acknowledge the support of Istituto Nazionale di Fisica Nucleare (INFN) iniziative specifiche MOONLIGHT2, QGSKY, and TEONGRAV. IDM acknowledges support from Grant IJCI2018-036198-I funded by MCIN/AEI/10.13039/501100011033 and, as appropriate, by “ESF Investing in your future” or by “European Union NextGenerationEU/PRTR”. IDM and RDM also acknowledge support from the grant PID2021-122938NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. DV also acknowledges the FCT project with ref. number PTDC/FIS-AST/0054/2021.