Joint constraints on cosmological parameters using future multi-band gravitational wave standard siren observations

The precise measurements of the cosmic microwave background (CMB) anisotropies initiated the era of precision cosmology [1, 2]. Nevertheless, with the improvement of measurement precisions of cosmological parameters, some tensions between the early- and late-universe observations have arised. In particular, the values of the Hubble constant inferred from the $ Planck $ CMB observation (based on the ΛCDM model) [3] and determined by the distance-ladder measurement (model-independent) [4] are shown to be in more than $ 5\sigma $ tension [4], which is now commonly believed as a severe crisis for cosmology [5, 6]. The Hubble tension has been widely discussed in the literature [5–30]. Thus far, there is no consensus on a valid extended cosmological model that can truly solve the Hubble tension. Therefore, some cosmological probes that can independently measure the Hubble constant need to be developed. The gravitational wave (GW) standard siren method is one of the most promising options.

Unlike the traditional electromagnetic (EM) observations, GW observations open a new window for exploring the expansion history of the universe. The GW waveform encodes the information of the luminosity distance, which is called a standard siren [31, 32]. Applying GW standard sirens in cosmology has recently been widely discussed in the literature [32–84]. If the redshift information of the GW source could be obtained by identifying the EM counterparts (we usually refer to this kind of GW standard sirens as bright sirens), the distance-redshift relation could be established for cosmological parameter estimations. While for the GW events without EM counterparts, the statistical analysis of the GW event associated with the galaxy catalog can also be applied to obtain the redshift information (we usually refer to this kind of GW standard sirens as dark sirens).

In fact, the frequency ranges of GW standard sirens are wide (from nanohertz to a few hundred hertz), corresponding to different GW sources. To detect GWs in different frequency bands, the pulsar timing arrays (PTAs), space-based GW detectors, and ground-based GW detectors are proposed.

The nanohertz GWs emitted by the supermassive black hole binaries (SMBHBs) could be detected by a PTA, a natural Galactic-scale detector of millisecond pulsars (MSPs). Although it is difficult to detect GWs from individual SMBHBs by the current PTA projects, e.g., the European Pulsar Timing Array [85], North American Nanohertz Observatory for Gravitational Waves [86], and Parkes Pulsar Timing Array (Australia) [87], it is expected that the individual SMBHBs could be detected by the SKA-era PTAs [88]. Yan et al. [89] proposed that the currently available SMBHB candidates with known redshifts could be detected by the future SKA-era PTAs, allowing SMBHBs to be treated as standard sirens to explore the cosmic expansion history. Wang et al. [62] forecasted the cosmological parameter estimation with bright and dark sirens of individual SMBHBs with the SKA-era PTAs.

The space-based GW detectors are proposed to detect GWs emitted from the massive black hole binaries (MBHBs) in the millihertz frequency band, e.g., Taiji [90–92], TianQin [93–95], and the Laser Interferometer Space Antenna [96–97]. The space-based GW detectors could detect high-redshift GW events (up to $ z\simeq 15-20 $), which are expected to provide high-redshift GW standard siren data. Some works show that the EM signals could be emitted in the process of MBHB mergers in both the radio and optical bands [98–105]. The applications of these bright sirens in cosmological parameter estimation have been forecasted in the literature [56, 57, 60, 97, 106–110].

The ground-based GW detectors could observe stellar-mass binaries in the frequency band of a few hundred hertz. The only multi-messenger observation event GW170817 from a binary neutron star (BNS) merger gave the first measurement of the Hubble constant using the standard siren method with a $ 14\% $ precision [111]. The measurement precision of the Hubble constant could reach $ 2\% $ using 50 similar GW standard sirens [112], showing the potential of the standard siren method in cosmological parameter estimation. While for the dark siren method, the latest constraint precision of the Hubble constant from the LIGO-Virgo-KAGRA observation is 19% [113] (for recent related works refer to, e.g., Refs. [114–116]). In the next decades, the third-generation (3G) ground-based GW detectors, the Cosmic Explorer (CE) in the U.S. [117] and the Einstein Telescope (ET) in Europe [118], will observe a large number of GW events in a wide range of redshift because their sensitivities are one order of magnitude improved over the current detectors [119].

In the next decades, it is expected that we could obtain the multi-band GW standard siren data. Owing to the fact that the numbers of detectable GWs and signal-to-noise ratios (SNRs) in different frequency bands are different, the joint future multi-band GW standard siren observations are expected to play an important role in cosmological parameter estimation.

In this work, first, the precision to which the cosmological parameters could be measured to by the joint constraints of future multi-band GW standard siren observations is determined. Second, the role of the multi-band GW standard sirens in breaking cosmological parameter degeneracies generated by the EM observations is discussed. Note that, in this work, we only focus on the GW bright standard siren observations. We will consider the future bright siren observations from the SKA-era PTAs, space-based GW detectors, and 3G ground-based GW detectors, which are in different frequency bands, and constrain the cosmological parameters relevant to dark energy and the Hubble constant issues using the mock data of joint multi-band GW standard sirens.

The remainder of this paper is organized as follows. In Sec. II.A, we introduce the method of simulating GW standard sirens from the SKA-era PTA. In Sec. II.B, we introduce the method of simulating GW standard sirens from Taiji. In Sec. II.C, we introduce the method of simulating GW standard sirens from CE. In Sec. III, we give the constraint results and make some relevant discussions. The conclusion is given in Sec. IV. We adopt the ΛCDM model as the fiducial model to generate the simulated GW standard siren data, with the cosmological parameters set to the constraint results obtained from $ Planck $ 2018 TT, TE, EE+lowE [3].

GW signals are detected in the timing residuals of MSPs by removing model-predicted times of arrival (ToAs) from the observational ToA data. The time residuals induced by a single GW source measured at time t on the Earth can be written as [89]

where $F_{+, \times}(\hat{\Omega}_{\rm s},\hat{\Omega}_{\rm p})$ are the antenna response functions [120], and $ \hat{\Omega}_{\rm s} $ and $ \hat{\Omega}_{\rm p} $ are the unit vectors pointing from the GW source ($ \alpha_{\rm s} $, $ \beta_{\rm s} $) and pulsar to the observer ($ \alpha_{\rm p} $, $ \beta_{\rm p} $), respectively. $ \Delta A_{+,\times}(t)=A_{+,\times}(t)-A_{+,\times}(t_{\rm p}) $ is the difference between the earth term $ A_{+,\times}(t) $ and the pulsar term $ A_{+,\times}(t_{\rm p}) $, where $ t_{\rm p} $ is the time at which GW passes the MSP [121]. The forms of $ A_{+,\times}(t) $ are related to the GW strain. Assuming SMBHBs inspiral in circular orbits, the GW strain $ h(t) $ can be written as

where $ \mathcal{M}_{\rm c} $ is the observed chirp mass, $ M=m_1+m_2 $ is the total mass of a binary system with the component masses $ m_1 $ and $ m_2 $, $ \eta=m_1 m_2/(m_1+m_2)^2 $ is the symmetric mass ratio, $ f_0=2f_{\rm orb} $ is the GW frequency at the time of our first observation, $ f_{\rm orb}=(2\pi T)^{-1} $ is the orbit frequency, and T is the orbital periods of SMBHB candidates taken from Refs. [122–129]. Here, we calculate $ f_{0} $ using the orbital periods of the 154 SMBHB candidates.

The SNR of the GW signal detected by a PTA is written as [89]

where $ N_{\rm p} $ is the number of MSPs, N is the total number of data points of each MSP, $ s_{i}(t_{n}) $ is the timing residual of i-th MSP at time $ t_{n} $, and $ \sigma_{t,i} $ is the root mean square (rms) timing residual of the i-th MSP. Here, we set the threshold of SNR to be 8.

We use the Fisher information matrix to estimate the measurement errors of $ d_{\rm L} $. For a PTA containing $ N_{\rm p} $ independent MSPs, the Fisher matrix is expressed as [89]

where $ \mathit{\boldsymbol{\theta}} $ denotes the free parameters to be estimated. Here, the Fisher matrix includes nine parameters, including eight GW source parameters ($ d_{\rm L} $, $ \mathcal{M}_{\rm c} $, $ \alpha_{\rm s} $, $ \beta_{\rm s} $, $\iota $, ψ, $ \phi_{0} $, $ f_{0} $) and the pulsar distance $ d_{\rm p} $. The error of the parameter $ \theta_a $ is calculated by $ \Delta\theta_a=\sqrt{(F^{-1})_{aa}} $, i.e., $ \sigma_{d_{\rm L}}^{\rm inst}=\Delta d_{\rm L}=\Delta\theta_1 $.

The measurement of $ d_{\rm L} $ is also affected by the weak lensing, and we adopt the form [130–132]

The total error of $ d_{\rm L} $ can be written as $\sigma_{d_{\rm L}}= \sqrt{(\sigma_{d_{\rm L}}^{\rm inst})^2+(\sigma_{d_{\rm L}}^{\rm lens})^2}$.

We analyze a catalog of 154 currently available SMBHB candidates for this work. Among these, 149 are obtained through periodic variations in their light curves, as described in previous studies [123–125]. The remaining candidates include Mrk 231 from Ref. [126], NGC 5548 from Ref. [127], OJ 287 from Ref. [122], SDSS J0159+0105 from Ref. [128], and Ark 120 from Ref. [129]. We obtain the right ascension, declination, redshift, and total mass information for the SMBHB candidates and fix other parameters. A previous work has shown that the polarization angle ψ and initial phase $ \phi_0 $ have no significant effects on the GW analysis [89]. Therefore, we follow Refs. [62, 89] and assume $ \psi = 0^{\circ} $ and $ \phi_0 = 0^{\circ} $ in our simulation. For the inclination angle $\iota $, we assume that all the GW events have an edge-on inclination angle, i.e., $ \iota=90^{\circ} $. While this assumption may not always hold true, we adopt it as a conservative analysis. In addition, note that we make an optimistic assumption that the mass ratios of the SMBHB candidates are assumed to be $ q=1 $, which is also adopted in Refs. [62, 89]. In fact, the mass ratio can impact the strength of the GW signal, with smaller values of q leading to weaker GW signals, larger measurement errors of $ d_{\rm L} $, and worse constraint results, as discussed in previous literature (e.g., Ref. [89]).

Here, we note that the ability of PTA GW observations is affected by many factors, for example, the number of MSPs $ N_{\rm p} $ and rms of time residual $ \sigma_t $. It is found that about 100 high-quality MSPs are sufficient for detecting nanohertz GWs from individual SMBHBs [62]. The current PTAs usually contain dozens of MSPs. In the future, we can expect SKA and FAST to observe hundreds of MSPs. In the present work, we simulate 200 pulsars to detect nanohertz GWs. For the rms of timing residual, it consists of white noise and red noise. Recent analysis shows that the total white noises of pulsars could approach 10–50 ns [133] in the SKA era; thus, we consider $ \sigma_{t} $ to be $ \sigma_{t}=20 $ ns. SMBHB candidates usually emit GWs in the frequency range of $ 10^{-7} $–$ 10^{-8} $ Hz. In this frequency, the red noise can be attenuated to a low noise level, and it does not affect the single GW detection, so it can be ignored. In addition, the stochastic gravitational wave background (SGWB) will also affect the detection of SMBHBs. However, recent studies have shown that SGWB is likely to be detected in about ten years. The SGWB can be regarded as red noise, which has a slight impact in the frequency of $ 10^{-7} $–$ 10^{-8} $ Hz. Following Ref. [62], we assume that the ToA data are obtained by monitoring the MSPs with the typical cadence of two weeks and the observation time is 10 years. Based on the simulation method introduced above, we simulate 35 bright sirens for the 10-year observation of SKA-era PTAs. Note that the number of detected SMBHBs in the SKA era is expected to be much larger than 154. A full analysis of the expected detection number of SMBHBs in the SKA era is left to a future work.

In this section, we focus on the GW signal from the inspiral of a non-spinning MBHB. The frequency domain GW waveform is written as

The effective luminosity distance, denoted as $ D_{\rm eff} $, is given by the formula $ D_{{\rm eff}}=d_{\rm L} [F^{2}_{+}(\frac{1+{\rm cos}^2 \iota}{2})^2+F^{2}_{\times} {\rm cos}^2 \iota]^{-1/2} $. Here, $ d_{\rm L} $ is the luminosity distance, and $ F_{+,\times} $ are the antenna response functions that depend on the location of the GW source (θ, ϕ) and the polarization angle ψ. The detailed expressions for $ F_{+,\times} $ can be found in Ref. [134]. The GW phase Ψ is written to the second Post-Newtonian order and is related to the coalescence time $ t_{\rm c} $ and coalescence phase $ \psi_{\rm c} $. The specific form of Ψ can be found in Ref. [134]. To describe the GW signal in Fourier space, the observation time t is replaced with $t(f)= t_{\rm c}-\frac{5}{256}\mathcal{M}^{-5/3}_{\rm c}(\pi f)^{-8/3}$ [135, 136].

The combined SNR for the detection network of N independent interferometers is given by

where $ \rho_i=\sqrt{\langle\tilde{h_i},\tilde{h_i}\rangle} $. The inner product is defined as

where $ f_{\rm lower}=10^{-4} $ Hz and $ f_{\rm upper}=c^{3}/6\sqrt{6}\pi GM_{\rm obs} $ with $ M_{\rm obs}=(m_{1}+m_{2})(1+z) $ [137]. Taiji's PSD is taken from Ref. [91]. We set the threshold of SNR to be 8 in the simulation.

For a network with N independent interferometers, the Fisher matrix can be written as

where $ \mathit{\boldsymbol{\theta}} $ denotes nine GW source parameters ($ d_{\rm L} $, $ \mathcal{M}_{\rm c} $, η, $ t_{\rm c} $, $ \psi_{\rm c} $, $\iota $, θ, ϕ, ψ). The instrumental error of the luminosity distance is $ \sigma_{d_{\rm L}}^{\rm inst}=\Delta d_{\rm L}=\sqrt{(F^{-1})_{11}} $. We use Eq. (6) to calculate the weak-lensing error. The error caused by the peculiar velocity of the GW source is given by [138]

where $ H(z) $ is the Hubble parameter and $ \sqrt{\langle v^{2}\rangle} $ is the peculiar velocity of the GW source. In this work, we set $ \sqrt{\langle v^{2}\rangle}=500 $ km/s, in agreement with the average values observed in galaxy catalogs [131]. In addition, we make the assumption that redshift measurements at $ z<2 $ are determined spectroscopically with negligible errors (Refs. [139, 140]). However, for the GW event with $ z>2 $ associated with photometric measurements, the redshift error is estimated as $ \Delta z\approx 0.03(1+z) $ [141, 142]. Following Refs. [56, 57, 106, 130, 131], we propagate the redshift error to the distance error by assuming our fiducial cosmology, i.e., $ \sigma^{\rm reds}_{d_{\rm L}}=\dfrac{\partial d_{\rm L} }{\partial z}\Delta z $. In fact, as shown in Fig. 1, for the GW events at $ z>2 $, the error from lensing is dominant. Therefore, the treatment has little effect on the cosmological parameter estimation. The total error of $ d_{\rm L} $ can be written as $\sigma_{d_{\rm L}}= \sqrt{(\sigma_{d_{\rm L}}^{\rm inst})^2+(\sigma_{d_{\rm L}}^{\rm lens})^2+(\sigma_{d_{\rm L}}^{\rm pv})^2+(\sigma_{d_{\rm L}}^{\rm reds})^2}$.

Figure 1.  (color online) Relative errors of luminosity distances from lensing, peculiar velocities, redshift measurements (only for GW standard sirens at $ z>2 $ from Taiji), and instrumental errors of PTA, Taiji, and CE.

As the origin of MBHs is currently unclear, there is uncertainty in predicting the event rate of MBHBs. Based on a semianalytical galaxy formation model, three population models of MBHBs, the pop III, Q3d, and Q3nod models are proposed. The three models have different mechanisms of seeding and delay. In fact, previous works [57, 130] have shown that the three population models can lead to different constraints on cosmological parameters due to the difference in the number of standard sirens based on the three models. The pop III model typically offers intermediate constraints among the three cases. Therefore, in the present work, we generate simulated standard siren data based on the pop III model. Some works show that the EM signals could be emitted in the process of MBHB mergers in both the radio and optical bands [98–105]. Recent works predicted the number of GW detected by space-based GW observatory where the EM counterparts could also be detected by SKA, ELT, and LSST [57, 108, 130, 143, 144]. Following Refs. [57, 108, 130, 143, 144], we select the GW events with $ \rm SNR>8 $ and the sky localization error of $ \Delta \Omega<10\; \rm deg^2 $ (corresponding to the field of view of LSST) as potential candidates for EM counterpart detections. A previous work has shown that the EM counterparts detected by LSST are also detectable for SKA+ELT [130]. Therefore, in the present work, we only consider the case of SKA+ELT. In this case, EM counterparts may first be detected by SKA in the radio band and the host galaxies are then identified through localization [130]. Then, the redshifts are determined spectroscopically or photometrically by the optical telescope ELT. We calculate the radio and optical luminosities of an EM counterpart, and it can be detected if its luminosities meet the thresholds of SKA and ELT, as discussed in, e.g., Refs. [57, 130, 144].

For each simulated standard siren event, the sky location (θ, ϕ), binary inclination $\iota $, coalescence phase $ \psi_{\rm c} $, and polarization angle ψ are evenly sampled in the ranges of $ \cos\theta \in [-1, 1] $, $ \phi \in [0,360^{\circ}] $, $ \cos\iota \in [-1, 1] $, $ \psi_{\rm c} \in $ [0, $ 360^\circ $], $ \psi \in $ [0, $ 360^\circ $], respectively. In this work, we assume $ t_{\rm c}=0 $ for simplicity. For the redshift and mass distributions of MBHBs, we use the numerical fitting formulas in Ref. [57] to fit the curves shown in Fig. 3 of Ref. [145]. According to the simulation method introduced above, we simulate 28 bright sirens for the five-year observation of Taiji based on the pop III model, which are shown in Fig. 2. The number is also basically consistent with that given in Ref. [57].

Figure 2.  (color online) Simulated GW standard siren data points observed by PTA, Taiji, and CE. The yellow data points represent the 35 standard sirens from the 10-year observation of PTA, red data points represent the 28 standard sirens from the 5-year observation of Taiji, and blue data points represent the 1000 standard sirens from the 10-year observation of CE.

In this work, we consider that all the GW standard sirens that can be detected by CE are the BNS mergers. The redshift distribution of BNS mergers adopts the form in Refs. [39, 45]. we adopt the GW waveform of the inspiralling non-spinning BNS system [45]. We use Eq. (8) to calculate the SNR of each GW event. For CE, $ f_{\rm lower}= $ 1 Hz is the lower cutoff frequency and $f_{\rm upper}=2/ $ $ (6^{3/2}2\pi M_{\rm obs})$ is the frequency at the last stable orbit with $ M_{\rm obs}=(m_1+m_2)(1+z) $. We adopt the PSD of 40 km-arm-length CE [146]. Here, we set the threshold of SNR to 8 in our simulation. In this work, we consider three measurement errors of $ d_{\rm L} $, including the instrumental error $ \sigma^{\rm inst}_{d_{\rm L}} $, weak-lensing error $ \sigma^{\rm lens}_{d_{\rm L}} $, and peculiar velocity error $ \sigma^{\rm pv}_{d_{\rm L}} $.

Recent forecasts show that the 3G ground-based GW detectors would detect $ \mathcal{O}(10^5) $ BNS mergers per year, but only approximately $ 0.1\% $ of them have the detectable EM counterparts [147]. Chen et al. recently showed that 910 GW standard sirens could be detected based on the 10-year observation of CE and Swift++ [75]. Therefore, in the forecast in the present work, we simulate 1000 GW standard sirens generated by BNS mergers based on the 10-year observation of CE.

For each simulated standard siren event, the masses of NSs ($ m_{1} $, $ m_{2} $) are randomly chosen in the ranges of [1, 2] $ M_{\odot} $. Without loss of generality, the merger time is chosen to $ t_{\rm c}=0 $ in our analysis. Here, we note that the inclination angle should be randomly chosen in the range of $ \cos\iota \in [-1, 1] $ when simulating isotropic GW sources. However, in this work, we assume that the redshifts of the GW events are determined by detecting SGRBs. Since SGRBs are strongly beamed, the detectable inclination angle is approximately $ \iota\leq 20^{\circ} $ [75, 131, 132, 138, 148, 149]. Therefore, in the present work, $ \iota $ is randomly chosen in the range of $ \iota\in [0,20^\circ] $. The simulation of other angles is the same as that described in the simulation of MBHB. Based on the above analysis, we simulate 1000 bright sirens for the 10-year observation of CE, as shown in Fig. 2.

In Fig. 1, we present the relative errors of luminosity distances due to lensing, peculiar velocities, redshift measurements (only for GW standard sirens at $ z>2 $ from Taiji), and instrumental error of PTA, Taiji, and CE. For PTA, both the instrumental error and lensing error dominate the error of $ d_{\rm L} $. However, for Taiji, the error of $ d_{\rm L} $ is dominated by lensing owing to the high SNRs of GW events from Taiji. For CE, the instrumental error dominates the error of $ d_{\rm L} $. We combine the errors for every standard siren and construct the Hubble diagram using the future multi-band standard siren observations. In Fig. 2, we show the simulated standard siren data from PTA, Taiji, and CE. We can see that the number of standard sirens from CE is the largest, followed by PTA and Taiji. However, because the SNRs of GW events observed by CE are lower than those of PTA and Taiji, CE has the largest errors of luminosity distances at the same redshifts. The luminosity distance errors of PTA and Taiji are almost the same at similar redshifts.

We adopt the Markov Chain Monte Carlo method [150] to maximize the likelihood $ \mathcal{L}\propto\exp(-\chi^2/2) $ and infer the posterior probability distributions. The $ \chi^2 $ function is defined as

where N is the number of standard siren data points.

In this section, we report the constraint results. We use the simulated GW standard siren data from PTA, Taiji, and CE to constrain the ΛCDM [$ w(z)=-1 $], wCDM [$ w(z)=\rm constant $], and $ w_0w_a $CDM [$ w(z)=w_0+w_a z/(1+z) $] models by performing the Markov-chain Monte Carlo analysis [150]. For the CMB data, we employ the "$ Planck $ distance priors" from the $ Planck $ 2018 observation [3, 151]. The $ 1\sigma $ and $ 2\sigma $ posterior distribution contours for the cosmological parameters of interest are shown in Figs. 3–6, and the $ 1\sigma $ errors for the marginalized parameter constraints are summarized in Table 1. We use $ \sigma(\xi) $ and $ \varepsilon(\xi) $ to represent the 1σ absolute and relative errors of the parameter $\xi $, with $ \varepsilon(\xi) $ defined as $ \varepsilon(\xi)=\sigma(\xi)/\xi $.

Figure 3.  (color online) Constraints on the ΛCDM model. Left panel: Two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the $\Omega_\rm {m}$–$ H_0 $ plane by using the PTA, Taiji, CE, and PTA+Taiji+CE data. Right panel: Two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the $\Omega_\rm{m}$–$ H_0 $ plane by using the CMB, PTA+Taiji+CE, and CMB+PTA+Taiji+CE data.

Figure 4.  (color online) Constraints on the wCDM model. Here, we show the two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the w–$\Omega_\rm{m}$ (left panel) and w–$ H_0 $ (right panel) planes using the PTA, Taiji, CE, and PTA+Taiji+CE data.

Figure 5.  (color online) Constraints on the wCDM model. Here, we show the two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the w–$\Omega_\rm{m}$ (left panel) and w–$ H_{0} $ (right panel) planes using the CMB, PTA+Taiji+CE, and CMB+PTA+Taiji+CE data.

Figure 6.  (color online) Constraints on the $ w_0w_a $CDM model. Left panel: Two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the $ w_0 $–$ w_a $ plane by using the PTA, Taiji, CE, and PTA+Taiji+CE data. Right panel: Two-dimensional marginalized contours ($ 68.3\% $ and $ 95.4\% $ confidence level) in the $ w_0 $–$ w_a $ plane by using the CMB, PTA+Taiji+CE, and CMB+PTA+Taiji+CE data.

Model	Error	CMB	PTA	Taiji	CE	PTA+Taiji+CE	CMB+PTA+Taiji+CE
ΛCDM	$\sigma(\Omega_\rm m)$	0.009	0.020	0.024	0.012	0.008	0.003
	$ \sigma(H_0) $	0.61	0.49	0.70	0.45	0.29	0.20
	$\varepsilon(\Omega_\rm m)$	$ 2.7\% $	$ 6.3\% $	$ 11.7\% $	$ 3.8\% $	$ 2.8\% $	$ 1.0\% $
	$ \varepsilon(H_0) $	$ 0.9\% $	$ 0.7\% $	$ 1.6\% $	$ 0.7\% $	$ 0.5\% $	$ 0.3\% $
wCDM	$\sigma(\Omega_\rm m)$	0.058	0.038	0.033	0.017	0.015	0.003
	$ \sigma(H_0) $	6.30	2.25	1.75	0.95	0.69	0.37
	$ \sigma(w) $	0.210	0.395	0.235	0.120	0.101	0.028
	$\varepsilon(\Omega_\rm m)$	$ 18.5\% $	$ 11.8\% $	$ 10.2\% $	$ 5.4\% $	$ 4.7\% $	$ 1.0\% $
	$ \varepsilon(H_0) $	$ 9.1\% $	$ 3.3\% $	$ 2.6\% $	$ 1.4\% $	$ 1.0\% $	$ 0.6\% $
	$ \varepsilon(w) $	$ 20.2\% $	$ 31.9\% $	$ 22.0\% $	$ 11.8\% $	$ 10.0\% $	$ 2.8\% $
$ w_{0}w_{a} $CDM	$\sigma(\Omega_\rm m)$	0.066	0.069	0.059	0.048	0.047	0.009
	$ \sigma(H_0) $	7.25	4.30	3.10	1.40	1.35	0.85
	$ \sigma(w_{0}) $	0.605	0.750	0.530	0.220	0.195	0.110
	$ \sigma(w_{a}) $	2.50	2.90	2.80	1.31	1.22	0.32
	$\varepsilon(\Omega_\rm m)$	$ 20.6\% $	$ 19.9\% $	$ 16.9\% $	$ 14.8\% $	$ 14.8\% $	$ 2.7\% $
	$ \varepsilon(H_0) $	$ 10.6\% $	$ 6.4\% $	$ 4.7\% $	$ 2.1\% $	$ 2.0\% $	$ 1.3\% $
	$ \varepsilon(w_{0}) $	$ 112.0\% $	$ 76.5\% $	$ 74.6\% $	$ 24.2\% $	$ 21.0\% $	$ 10.9\% $

Table 1.  Absolute errors (1σ) and relative errors of the cosmological parameters in the ΛCDM, wCDM, and $ w_{0}w_{a} $CDM models using the CMB, PTA, Taiji, CE, PTA+Taiji+CE, and CMB+PTA+Taiji+CE data. Here, $ H_0 $ is in units of km s$ ^{-1} $ Mpc$ ^{-1} $.

We first focus on the constraint results for the ΛCDM model. In the left panel of Fig. 3, we show the constraint results in the $\Omega_\rm {m}$–$ H_0 $ plane by using the simulated PTA, Taiji, CE, and PTA+Taiji+CE data. As can be seen, CE contributes the most to the PTA+Taiji+CE results, followed by PTA and Taiji. This is because the number of simulated standard sirens from CE is much higer than that for PTA and Taiji. Although the measurement errors of $ d_{\rm L} $ for CE are large, the constraints on cosmological parameters are reduced statistically. Compared with Taiji, PTA has more lower-redshift data points ($ z<2 $), so PTA can better constrain the Hubble constant. Moreover, owing to the different redshift intervals of the simulated PTA, Taiji, and CE data, their parameter degeneracy orientations are slightly different, so their combination could break the cosmological parameter degeneracies. The combination of PTA, Taiji, and CE gives $\sigma(\Omega_{\rm {m}})=0.008$ and $ \sigma(H_0)=0.29 $ km s$ ^{-1} $ Mpc$ ^{-1} $, which are $ 33.3\% $ [$(0.012- 0.008)/0.012$] and $ 35.6\% $ [$ (0.45-0.29)/0.45 $] better than those of CE. In the right panel of Fig. 3, we can see that the contours of CMB and PTA+Taiji+CE show different orientations; thus, their combination could break the cosmological parameter degeneracies. The prime cause is that GW could measure $ H_0 $ better, so it could lead to a different degeneracy direction compared with CMB. The combination of CMB and PTA+Taiji+CE gives $\sigma(\Omega_{\rm{m}})=0.003$ and $ \sigma(H_0)=0.20 $ km s$ ^{-1} $ Mpc$ ^{-1} $, and the constraint precisions of $\Omega_\rm {m}$ and $ H_0 $ are $ 1.0\% $ and $ 0.3\% $, respectively, meeting the standard of precision cosmology. In general, the joint PTA+Taiji+CE data could tightly constrain the Hubble constant, and if combined with CMB, the measurement precisions of cosmological parameters could be greatly improved (better than or at least equal to $ 1\% $), because of the breaking of the parameter degeneracies.

In Figs. 4 and 5, we show the constraint results in the w–$\Omega_\rm{m}$ and w–$ H_0 $ planes for the wCDM model. As can be seen from Fig. 4, CE also contributes the most to the PTA+Taiji+CE data. However, the ability of Taiji to constrain the cosmological parameters in the wCDM model is better than that of PTA. This is because Taiji has more high-redshift standard sirens than PTA, so Taiji can better constrain the EoS parameter of dark energy w. Meanwhile, the combination of PTA, Taiji, and CE gives $ \sigma(w)=0.101 $, which is $ 15.8\% $ better than the constraint result by CE. In Fig. 5, we see that the parameter degeneracy orientations of CMB and PTA+Taiji+CE are almost orthogonal; therefore, their combination could not only the break cosmological parameter degeneracies but also tremendously improve the cosmological parameter constraints. The addition of the PTA+Taiji+CE data could reduce the $ 1\sigma $ absolute error of w by $ 86.7\% $, compared with CMB. Moreover, the combination of CMB and PTA+Taiji+CE gives $ \sigma(w)=0.028 $, which is comparable with the latest constraint result by the CMB+BAO+SN data [152].

In Fig. 6, we show the case for the $ w_0w_a $CDM model in the $ w_0 $–$ w_a $ plane. As can be seen from the left panel of Fig. 6, the constraint results are the same as those in the wCDM model, i.e., CE contributes the most, followed by Taiji and PTA. The joint PTA+Taiji+CE data could give $ \sigma(w_0)=0.195 $ and $ \sigma(w_0)=1.22 $, which are both better than the constraint results by the CMB data. Furthermore, in the right panel of Fig. 6, we see that the combination of CMB and PTA+Taiji+CE could also break the parameter degeneracies, thus significantly improving the cosmological parameter constraints. The joint CMB+PTA+Taiji+CE data give $ \sigma(w_0)=0.110 $ and $ \sigma(w_a)=0.32 $, which are $ 81.8\% $ [$ (0.605-0.110)/0.605 $] and $ 87.2\% $ [$ (2.50-0.32)/2.50 $] better than the results of CMB.

Our results show that the joint future multi-band GW standard siren observations would play a crucial role in cosmological parameter estimation. CE contributes the most to the PTA+Taiji+CE results since the number of standard sirens detected by CE is considerably higher than those of PTA and Taiji. PTA has more lower-redshift ($ z<2 $) data. Taiji has more higher-redshift data. Hence, PTA offers better constraints in the ΛCDM model, while in the dynamical dark energy models, Taiji offers better constraints. The joint multi-band GW standard siren data show great potential in constraining the ΛCDM model. Moreover, their parameter degeneracy orientations are slightly different, and thus, their combination could break the cosmological parameter degeneracies. However, the joint constraints do not perform well in the wCDM and $ w_0w_a $CDM models. Fortunately, the joint PTA+Taiji+CE data have different parameter degeneracy orientations from CMB; theorefore, their combination could effectively break the parameter degeneracies and greatly improve the constraint precisions of the cosmological parameters. It can be concluded that the future multi-band GW observations can be expected to aid in the precise measurement of the cosmological parameters and to help solve important cosmological problems. Here, we emphasize that the successful application of the standard siren method heavily depends on the accuracy of the GW data [153, 154]. In particular, the systematic errors caused by calibrations in the data processing should be carefully avoided. In this work, the systematic errors from calibrations in the standard siren data are not considered (for the impact of calibration uncertainties on the Hubble constant measurements, see e.g. Ref. [155]).

In this work, we explore the potential of the joint constraints on the cosmological parameters using future multi-band GW standard siren observations. We simulated the multi-band standard siren data based on the 10-year observation of CE, the 5-year observation of Taiji, and the 10-year observation of the SKA-era PTA, and we used mock data to constrain three typical cosmological models, i.e., the ΛCDM, wCDM, and $ w_0w_a $CDM models.

We find that the joint PTA+Taiji+CE data could provide tight constraints on the Hubble constant, with the constraint precision being $ 0.5\% $ in the ΛCDM model. However, the joint data do not perform well in constraining the EoS parameters of dark energy. Fortunately, CMB and PTA+Taiji+CE show different parameter degeneracy orientations; thus, their combination could effectively break the parameter degeneracies and improve constraint precisions of cosmological parameters. In the ΛCDM model, the constraint precisions of $\Omega_\rm{m}$ and $ H_0 $ using the CMB+PTA+Taiji+CE data are better than or at least equal to $ 1\% $. While in the wCDM model, CMB+ PTA+Taiji+CE offers $ \sigma(w)=0.028 $, which is comparable with the latest constraint result by the CMB+BAO+SN data. Compared with CMB, the combination of CMB and PTA+Taiji+CE could improve the constraint on w by $ 86.7\% $. In the $ w_0w_a $CDM model, the CMB+PTA+Taiji+ CE data offer $ \sigma(w_0)=0.110 $ and $ \sigma(w_a)=0.32 $, which are $ 81.8\% $ and $ 87.2\% $ better than the results by CMB.

Hence, we can conclude the following: (i) the joint future multi-band GW standard sirens could precisely measure the Hubble constant, but are not good at measuring dark energy; (ii) the joint PTA+Taiji+CE data could effectively break the cosmological parameter degeneracies generated by the CMB data, especially in the dynamical dark energy models. The future multi-band GW observations can be expected to be used to probe the nature of dark energy and measure the Hubble constant.

We thank Yong Yuan, Tao Han, and Peng-Ju Wu for helpful discussions.