This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Prospects for measuring the Hubble constant and dark energy using gravitational-wave dark sirens with neutron star tidal deformation

Shang-Jie Jin    Tian-Nuo Li    Jing-Fei Zhang    and Xin Zhang11footnotetext: Corresponding author.
Abstract

Using the measurements of tidal deformation in the binary neutron star (BNS) coalescences can obtain the information of redshifts of gravitational wave (GW) sources, and thus actually the cosmic expansion history can be investigated using solely such GW dark sirens. To do this, the key is to get a large number of accurate GW data, which can be achieved with the third-generation (3G) GW detectors. Here we wish to offer an answer to the question of whether the Hubble constant and the equation of state (EoS) of dark energy can be precisely measured using solely GW dark sirens. We find that in the era of 3G GW detectors 𝒪​(105−106){\cal O}(10^{5}-10^{6}) dark siren data (with the NS tidal measurements) could be obtained in three-year observation if the EoS of NS is perfectly known, and thus using only dark sirens can actually achieve the precision cosmology. Based on a network of 3G detectors, we obtain the constraint precisions of 0.15%0.15\% and 0.95%0.95\% for the Hubble constant H0H_{0} and the constant EoS of dark energy ww, respectively; for a two-parameter EoS parametrization of dark energy, the precision of w0w_{0} is 2.04%2.04\% and the error of waw_{a} is 0.13. We conclude that 3G GW detectors would lead to breakthroughs in solving the Hubble tension and revealing the nature of dark energy provided that the EoS of NS is perfectly known.

1 Introduction

Allan Sandage pointed out in as early as 1970 that the two key numbers that should be precisely measured in cosmology are the Hubble constant (H0H_{0}) and the deceleration parameter (q0q_{0}) [1]. This is undoubtedly a deep insight. The measurement of q0q_{0} using type Ia supernovae has led to the discovery of the accelerating expansion of the universe [2, 3], which is usually explained through assuming an exotic component with a negative pressure that is currently dominating the evolution of the universe, dubbed dark energy. However, the nature of dark energy is still an enigma, which might be solved if the equation of state (EoS) of dark energy is precisely measured. Therefore, in the current cosmology, an important mission is to precisely measure the Hubble constant and the EoS of dark energy.

The Hubble constant is a direct measure of the current expansion rate of the universe and it is a key cosmological parameter because all the cosmic properties relevant to the absolute scale, such as the age and size of the universe as well as some physical processes like the growth of cosmic structure and the nucleosynthesis of the light isotopes, are determined by it. Lots of efforts have been made to measure the Hubble constant in the past near one century, and currently the precisions of the measurements have reached (for indirect measurement) or nearly reached (for direct measurement) 1%1\%. However, the values of the Hubble constant inferred from the observations of the P​l​a​n​c​kPlanck cosmic microwave background (CMB) anisotropies (a 0.8%0.8\% measurement, assuming a flat Λ\LambdaCDM cosmology) [4] and determined from the Cepheid-supernova distance ladder (a 1.4%1.4\% measurement) [5] are highly inconsistent, with the tension reaching 5σ\sigma significance [5], which is the so-called “Hubble tension”, greatly challenging the current cosmology. Recently, the Hubble tension has been widely discussed in the literature (see, e.g., refs. [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]).

For the measurement of the EoS of dark energy (assuming a constant ww), the CMB+BAO+SN data have given the tightest constraint to date, with the constraint precision being about 3%3\% [32]. However, such a measurement precision is still far away from revealing the fundamental nature of dark energy.

The gravitational wave (GW) standard siren observation opened a new window into studying the expansion history of the universe, which recently has been widely discussed [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75] (see e.g., ref. [76] for a recent review). The luminosity distance of a GW source could be directly obtained through the analysis of the GW waveform without the external calibration, which is called a standard siren [77, 33]. Furthermore, some GW sources, e.g., binary neutron star (BNS) mergers, could produce rich electromagnetic (EM) signals, such as kilonovae, short γ\gamma-ray bursts and their afterglows, referred to as EM counterparts of GWs [78, 79]. GWs together with their EM counterparts, if can both be observed, realize the multi-messenger observations, which are usually referred to as bright sirens, and could enable us to establish the absolute luminosity distance–redshift relation to constrain H0H_{0} and w​(z)w(z). In fact, the only multi-messenger observation, GW170817 from the first detected BNS merger, has given an independent constraint on the Hubble constant, at a 14%14\% precision [80]. Limited by the sensitivities of the current GW detectors (LIGO and Virgo), only two BNS mergers are detected during the first, second, and third observing runs [81, 82, 83] and only a few BNS mergers will be detected during the fourth observing run [84]. The Hubble constant would be measured to a 2%2\% precision using about 50 GW170817-like standard sirens [85], which may help resolve the Hubble tension. However, in order to effectively measure the EoS of dark energy, considering only the second-generation GW detectors is definitely not hopeful, and one would has to resort to the next-generation detectors.

In the next decades, the third-generation (3G) ground-based GW detectors, the Cosmic Explorer (CE) [86] and the Einstein Telescope (ET) [87], will have a powerful detection sensitivity, one order of magnitude improved over the current GW detectors [88]. A series of forecasts [89, 90, 91] show that the 3G GW detectors will detect 𝒪​(105)\mathcal{O}(10^{5}) BNS mergers per year. However, even so, only a small part of them (around 0.1%0.1\%) have the detectable EM counterparts [91], and thus the application of bright sirens in cosmology has strong limitations. Furthermore, it has been shown in some works [76, 52, 92, 93] that solely using GW bright sirens cannot tightly constrain the EoS of dark energy and the main role they play in cosmology is to help break the cosmological parameter degeneracies of other traditional observations.

In 2012, Messenger and Read [94] proposed that the distance–redshift relation can be established by solely using the GW observations based on the 3G GW detectors (note here that the GW standard sirens without EM counterparts are called “dark sirens”). This is because the tidal deformations of NSs in the BNS systems could be measured using the 3G GW detectors. The tidal deformations depend on the EoS of NS and could provide an additional contribution to the phase evolution of the GW waveform. The additional phase is related to the intrinsic mass and thus it can be used to break the degeneracy between mass and redshift, which leads to the redshift of BNS being measured.

Recently, using the GW dark sirens (tidal measurements of BNS) to measure the cosmological parameters has been discussed in the literature [95, 96, 97, 98]. Del Pozzo et al. [95] forecasted the ability of GW dark sirens from ET to measure the cosmological parameters. Wang et al. [96] forecasted the constraints on w0w_{0} and waw_{a} (with other cosmological parameters fixed) using GW dark sirens from the 3G ground-based GW network. Chatterjee et al. [97] forecasted the measurements of the Hubble constant using GW dark sirens from O5, Voyager, and CE. Ghosh et al. [98] forecasted the measurements of the Hubble constant using GW dark sirens from CE.

Previous works focused on the ability to measure the Hubble constant using the mock dark siren data of a single 3G GW observatory. Actually, in the next decades, the 3G GW detectors are expected to form a ground-based detector network, thereby improving the ability to detect BNS mergers. In this work, we wish to answer the question of whether the Hubble constant and the EoS of dark energy can both be precisely measured using only the GW dark sirens detected by the 3G ground-based GW detector network. We make cosmological analysis using four cases of 3G GW observations, single ET, single CE, CE-CE network (one CE in the US and another one in Australia, abbreviated as 2CE hereafter), and ET-CE-CE network (abbreviated as ET2CE hereafter). We note that GW dark sirens using the cross-correlation of binary black hole mergers and galaxy catalogs [99, 100, 101, 102, 61, 103, 104, 105, 55, 58] are not discussed here. For cosmological models, we consider the wwCDM and w0​waw_{0}w_{a}CDM models to make cosmological analysis. We employ Λ\LambdaCDM as the fiducial model to generate the mock GW dark siren data, with the fiducial values of cosmological parameters being set to the constraint results of P​l​a​n​c​kPlanck 2018 TT,TE,EE+lowE [4].

This work is organized as follows. In section 2.1, we briefly introduce the merger rate of BNS adopted in this work. In section 2.2, we introduce the method of calculating the detection number of BNS merger events. In section 2.3, we introduce the method of simulating GW dark sirens. In section 2.4, we briefly introduce the method of constraining cosmological parameters. In section 3, we show the constraint results. In section 4, we compare the constraint results of this work with those of other future observations. The conclusion is given in section 5.

2 Methodology

2.1 BNS merger rate

The merger rate as a function of redshift in the observer frame is given by

Ro​(z)=Rm​(z)1+z​d​V​(z)d​z,R_{\rm o}(z)=\frac{R_{\rm m}(z)}{1+z}\frac{dV(z)}{dz}, (2.1)

where the (1+z)(1+z) term arises from converting the source-frame time to the observer-frame time and d​V/d​zdV/dz is the comoving volume element. RmR_{\rm m} is the merger rate in the source frame, which is related to the cosmic star formation rate,

Rm​(zm)=∫zm∞𝑑zf​d​tfd​zf​Rsf​(zf)​P​(td),R_{\rm m}\left(z_{\rm m}\right)=\int_{z_{\rm m}}^{\infty}dz_{\rm f}\frac{dt_{\rm f}}{dz_{\rm f}}R_{\rm sf}\left(z_{\rm f}\right)P\left(t_{\rm d}\right), (2.2)

where d​t/d​z=−[(1+z)​E​(z)​H0]−1{dt}/{dz}=-[(1+z)E(z)H_{0}]^{-1}. The time-delay distribution P​(td)P(t_{\rm d}) is the probability density that the binary system formed at time tft_{\rm f} (or redshift zfz_{\rm f}) and merged at time tmt_{\rm m} (or redshift zmz_{\rm m}) with td=tm−tft_{\rm d}=t_{\rm m}-t_{\rm f}. RsfR_{\rm sf} is the cosmic star formation rate. We assume the Madau-Dickinson star formation rate [106], with an exponential time delay (an e-fold time of 100 Myr) [107]. Meanwhile, we adopt the local comoving merger rate to be 320 Gpc−3​yr−1\rm Gpc^{-3}~{}yr^{-1} that is the estimated median rate based on the O3a observation [108].

2.2 The detection number of BNS merger events

Then we calculate the detection number of BNS merger events. Here we adopt the signal-to-noise ratio (SNR) threshold to be 8 for a detection. The combined SNR for the detection network of NN detectors is ρ=∑a=1N(ρa)2\rho=\sqrt{\sum\limits_{a=1}^{N}(\rho_{a})^{2}}, where ρa\rho_{a} is the SNR of the aath detector and is given by

ρa2=56​(G​ℳchirp)5/3c3​π4/3​dL2​(z)​∫flowerfupper𝑑f​ℱa2​f−7/3Sa​(f),\rho_{a}^{2}=\frac{5}{6}\frac{\left(G\mathcal{M}_{\rm chirp}\right)^{5/3}}{c^{3}\pi^{4/3}d_{\rm L}^{2}(z)}\int_{f_{\rm lower}}^{f_{\rm upper}}df\frac{\mathcal{F}_{a}^{2}f^{-7/3}}{S_{a}(f)}, (2.3)

where ℳchirp=(1+z)​η3/5​M\mathcal{M}_{\rm chirp}=(1+z)\eta^{3/5}M is the observed chirp mass, M=m1+m2M=m_{1}+m_{2} is the total mass of a binary system with the component masses m1m_{1} and m2m_{2}, η=m1​m2/(m1+m2)2\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2} is the symmetric mass ratio, dLd_{\rm L} is the luminosity distance to the GW source, GG is the gravitational constant, cc is the speed of light, flowerf_{\rm lower} is the lower cutoff frequency (flower=1f_{\rm lower}=1 Hz for ET and flower=5f_{\rm lower}=5 Hz for CE), fupper=2/(63/2​2​π​Mobs)f_{\rm upper}=2/(6^{3/2}2\pi M_{\rm obs}) is the frequency at the last stable orbit with Mobs=(m1+m2)​(1+z)M_{\rm obs}=(m_{1}+m_{2})(1+z) [36], and Sa​(f)S_{a}(f) represents the one-side noise power spectral density (PSD) of the aa-th GW detector. We wish to note that in this work we consider the waveform in the inspiralling stage for a non-spinning BNS system. Hence, the upper cutoff frequency is chosen as the frequency at the last stable orbit. If the inspiral-merger-ringdown waveform is considered, a higher SNR is expected to provide better measurements on dLd_{\rm L}, thus providing better constraints on cosmological parameters. We adopt the PSD of ET from ref. [109] and that of CE from ref. [110] (we consider one 40 km-arm-length CE in the US and one 20 km-arm-length CE in Australia and the noise curves are shown in figure 1). ℱa=(1+cos2⁡ι)24​F+,a2+cos2⁡ι​F×,a2\mathcal{F}_{a}=\sqrt{\frac{\left(1+\cos^{2}\iota\right)^{2}}{4}F_{+,a}^{2}+\cos^{2}\iota F_{\times,a}^{2}} characterizes the detector response, where ι\iota is the inclination angle between the binary’s orbital angular momentum and the line of sight, F+,aF_{+,a} and F×,aF_{\times,a} are the aath detector’s antenna response functions of two GW polarizations, ++ and ×\times. Note that we adopt the NS’s mass distribution from the analysis of the latest O3 observation (figure 7 of ref. [111]).

Refer to caption
Figure 1: Sensitivity curves of ET [109] and CE [110] used in this work. Note that we consider one 40 km-arm-length CE in the US and one 20 km-arm-length CE in Australia.

The detector response functions are related to the source location, the polarization angle, the detector’s location (latitude φ\varphi, longitude λ\lambda, the angle between the interferometer arms ζ\zeta, and the orientation angle γ\gamma of the detector’s arms measured counter-clockwise from East of the Earth to the bisector of the interferometer arms). In table 1, we list the coordinates of the GW observatories considered in this work. Some other 3G GW detectors’ location candidates could be found in ref. [112]. We note that the 3G GW detectors have much wider frequency range than the current ground-based GW detectors. The time of BNS merger is a function of f−8/3f^{-8/3}. Thus, the BNS signals can be found in the 3G GW band for hours, even for several days, and in the 2G GW band for tens of minutes [113]. Therefore, the time dependence in the detector response functions is considered in this work. The forms of F+F_{+} and F×F_{\times} are related to these parameters, as detailedly described in refs. [114, 115, 116].

Table 1: The coordinates of the GW observatories [117, 118, 119, 120] considered in this work.
GW detector φ​[deg]\varphi\ [\mathrm{deg}] λ​[deg]\lambda\ [\mathrm{deg}] γ​[deg]\gamma\ [\mathrm{deg}] ζ​[deg]\zeta\ [\mathrm{deg}]
Cosmic Explorer, USA 43.82743.827 −112.825-112.825 45.00045.000 90.000
Einstein Telescope, Europe 40.44340.443 9.4579.457 0.0000.000 60.000
Cosmic Explorer, Australia −34.000-34.000 145.000145.000 90.00090.000 90.000
Refer to caption
Figure 2: The expected detection rate as a function of redshift for ET, CE, 2CE, and ET2CE.

According to our calculation, we obtain the numbers of the GW events detected by ET (545763), CE (1017618), the 2CE network (1317183), and the ET2CE network (1815501), assuming a three-year observation. We see that the number of GW events detected by CE is more than twice of ET. Compared to the single CE observatory, the number of GW events detected by the 2CE network almost doubles. Meanwhile, the number of GW events detected by the ET2CE network is slightly more than that of the 2CE network. In figure 2, we show the expected detection rate as a function of redshift for ET, CE, 2CE, and ET2CE.

2.3 GW dark sirens

Under the stationary phase approximation [121], the Fourier transform of the time-domain waveform for a detector network (with NdN_{d} detectors) is given by [113]

𝒉~​(f)=e−i​𝚽​𝒉^​(f),\tilde{\bm{h}}(f)=e^{-i\bf\Phi}\hat{\bm{h}}(f), (2.4)

with 𝒉^\hat{\bm{h}} given by

𝒉^​(f)=[h1​(f)S1​(f),h2​(f)S2​(f),⋯,hNd​(f)SNd​(f)],{\hat{\bm{h}}(f)=\left[\frac{{h}_{1}(f)}{\sqrt{S_{1}(f)}},\frac{{h}_{2}(f)}{\sqrt{S_{2}(f)}},\cdots,\frac{{h}_{N_{d}}(f)}{\sqrt{S_{N_{d}}(f)}}\right],} (2.5)

where 𝚽\bf\Phi is the Nd×NdN_{d}\times N_{d} diagonal matrix with Φa​b=2​π​f​δa​b​(n⋅ra)\Phi_{ab}=2\pi f\delta_{ab}(\textbf{n}\cdot\textbf{r}_{a}), n is the propagation direction of a GW, and ra\textbf{r}_{a} is the location of the aath detector. Here Sa​(f)S_{a}(f) is PSD of the aa-th detector. We consider the waveform in the inspiralling stage for a non-spinning BNS system. We adopt the restricted Post-Newtonian (PN) approximation and calculate the waveform to the 3.5PN order [122, 123],

ha​(f)=𝒜a\displaystyle h_{a}(f)=\mathcal{A}_{a} f−7/6​exp⁡{i​(2​π​f​tc−π/4−2​ψc+2​Ψ​(f/2)−φ(2,0),a+Φtidal​(f))},\displaystyle f^{-7/6}\exp\{i\big{(}2\pi ft_{\rm c}-\pi/4-2\psi_{\rm c}+2\Psi(f/2)-\varphi_{(2,0),a}+\Phi_{\text{tidal}}(f)\big{)}\}, (2.6)

where the Fourier amplitude 𝒜a\mathcal{A}_{a} is given by

𝒜a=\displaystyle\mathcal{A}_{a}= 1dL​(F+,a​(1+cos2⁡ι))2+(2​F×,a​cos⁡ι)2×5​π/96​π−7/6​ℳchirp5/6.\displaystyle~{}~{}\frac{1}{d_{\rm L}}\sqrt{\big{(}F_{+,a}(1+\cos^{2}\iota)\big{)}^{2}+(2F_{\times,a}\cos\iota)^{2}}\times\sqrt{5\pi/96}\pi^{-7/6}\mathcal{M}_{\rm chirp}^{5/6}. (2.7)

Here ψc\psi_{\rm c} is the coalescence phase and the detailed forms of Ψ​(f)\Psi(f) and φ(2,0),a\varphi_{(2,0),a} can be found in refs. [122, 123]. Note that we employ the stationary phase approximation to obtain the frequency-domain expression of GW waveform, in which the time tt in F+F_{+}, F×F_{\times}, and Φi​j\Phi_{ij} is replaced by tf=tc−(5/256)​ℳchirp−5/3​(π​f)−8/3t_{\rm f}=t_{\rm c}-(5/256)\mathcal{M}_{\rm chirp}^{-5/3}(\pi f)^{-8/3} [124, 113], with tct_{\rm c} the coalescence time.

The tidal contribution to the GW phase from a BNS system is given by [94]

Φtidal​(f)=∑j=123​λj128​η​M5​[−24χj​(1+11​ηχj)​(π​M​f)5/3−528​χj​(3179−919​χj−2286​χj2+260​χj3)​(π​M​f)7/3],\Phi_{\text{tidal}}(f)=\sum_{j=1}^{2}\frac{3\lambda_{j}}{128\eta M^{5}}\bigg{[}-\frac{24}{\chi_{j}}\big{(}1+\frac{11\eta}{\chi_{j}}\big{)}(\pi Mf)^{5/3}-\frac{5}{28\chi_{j}}\big{(}3179-919\chi_{j}-2286\chi_{j}^{2}+260\chi_{j}^{3}\big{)}(\pi Mf)^{7/3}\bigg{]}, (2.8)

where the sum is comprised of the components of the binary system, χj=mj/(m1+m2)\chi_{j}=m_{j}/(m_{1}+m_{2}), and λj\lambda_{j} is the tidal deformation, which is related to the EoS of NS and characterizes changes of the quadruple of NS. It is comparable in magnitude with the 3.5PN phasing term for NSs [94]. Due to the fact that EoS of NS is still unknown, we choose SLy [125] as the fiducial EoS of NS that is also consistent with the current observation [126]. In fact, the consideration of different EoSs of NS may lead to different constraint results. In ref. [94], Messenger and Read found that for the three typical NS’s EoSs, the redshift measurement errors obtained by using the tidal measurement method can differ by several times (a stiffer EoS gives better redshift measurement). In ref. [96], Wang et al. used different EoSs of NS and focused on constraining dark energy by fixing other cosmological parameters (a stiffer EoS of NS gives better constraints). Hence, in this work, in order to make the cosmological analysis representative, we choose a medium EoS of NS, namely SLy.

In this work, following ref. [96], in order to parameterize the EoS of NS, we express the λ\lambda–mm relation as a linear function of the NS mass in the range of 1.2​M⊙−2​M⊙1.2~{}M_{\odot}-2~{}M_{\odot} as

λ=B​m+C,\lambda=Bm+C, (2.9)

where BB and CC are tidal effect parameters. For the fitting of SLy, we obtain B=−2.14B=-2.14 and C=4.64C=4.64, where the goodness-of-fit value is 0.999. Therefore, the influence of the fitting errors on the results is almost negligible.

Then we can use the Fisher information matrix to calculate the measurement errors of the GW parameters. For a set of parameters 𝜽\bm{\theta}, the Fisher information matrix of a network including NdN_{d} independent GW detectors is given by

Fi​j=(∂𝒉~∂θi|∂𝒉~∂θj),{F}_{ij}=\left(\frac{\partial\tilde{\bm{h}}}{\partial\theta_{i}}\Bigg{|}\frac{\partial\tilde{\bm{h}}}{\partial\theta_{j}}\right), (2.10)

where θi\theta_{i} denotes twelve parameters (dLd_{\rm L}, zz, ℳchirp\mathcal{M}_{\rm chirp}, η\eta, tct_{\rm c}, ψc\psi_{\rm c}, ι\iota, δ\delta, α\alpha, ψ\psi, BB, CC) for a given BNS system. The inner product is defined as

(ζ|ξ)=2​∫flower fupper{ζ​(f)​ξ∗​(f)+ζ∗​(f)​ξ​(f)}​𝑑f,\left(\zeta|\xi\right)=2\int_{f_{\text{lower }}}^{f_{\mathrm{upper}}}\{{\zeta}(f){\xi}^{*}(f)+{\zeta}^{*}(f){\xi}(f)\}df, (2.11)

where ∗* represents complex conjugate. Here we numerically calculate the partial derivatives ∂h~​(f)/∂θi\partial\tilde{h}(f)/\partial\theta_{i} by [h~​(f,θi+d​θi)−h~​(f,θi)]/d​θi[\tilde{h}(f,\theta_{i}+d\theta_{i})-\tilde{h}(f,\theta_{i})]/d\theta_{i}, with d​θi=10−nd\theta_{i}=10^{-n}. In order to verify the robustness of our code, we compare our Fisher matrix with that obtained from the state-of-the-art code 𝙶𝚆𝙵𝙰𝚂𝚃\tt GWFAST [127]. We find that the results are consistent after optimizing the value of nn to ensure the convergence of the derivative. The Fisher information matrix can then be used to calculate the measurement errors of the GW parameters

Δ​θi=(F−1)i​i.\Delta\theta_{i}=\sqrt{(F^{-1})_{ii}}. (2.12)

The instrumental error is σdLinst=Δ​dL=(F−1)11\sigma_{\rm d_{\rm L}}^{\rm inst}=\Delta d_{\rm L}=\sqrt{(F^{-1})_{11}}. The measurement of dLd_{\rm L} is also affected by the weak lensing and we adopt the form [128, 129, 130]

σdLlens​(z)=[1−0.3π/2​arctan⁡(z/0.073)]×dL​(z)×0.066​[1−(1+z)−0.250.25]1.8.\sigma_{d_{\rm L}}^{\rm lens}(z)=\left[1-\frac{0.3}{\pi/2}\arctan\left(z/0.073\right)\right]\times d_{\rm L}(z)\times 0.066\bigg{[}\frac{1-(1+z)^{-0.25}}{0.25}\bigg{]}^{1.8}. (2.13)

The total error of dLd_{\rm L} is given by

σdL=(σdLinst)2+(σdLlens)2.\displaystyle\sigma_{d_{\rm L}}=\sqrt{(\sigma_{d_{\rm L}}^{\rm inst})^{2}+(\sigma_{d_{\rm L}}^{\rm lens})^{2}}. (2.14)
Refer to caption
Figure 3: Percentage of Δ​z/z\Delta z/z distribution of GW events detected by ET, CE, 2CE, and ET2CE. In the upper four panels, we show the percentages of Δ​z/z\Delta z/z distributions for ET, CE, 2CE, and ET2CE in blue, orange, red, and green lines, respectively. In the lower panel, we use the interpolation method to fit the upper four panels to show the distribution of Δ​z/z\Delta z/z intuitively.

Throughout this work, we simulate the BNS mergers with random binary orientations and sky directions. The sky direction, binary inclination, polarization angle, and the coalescence phase are randomly chosen in the ranges of cos​δ∈[−1,1]{\rm cos}\,\delta\in[-1,1], α∈[0,360∘]\alpha\in[0,360^{\circ}], cos​ι∈[−1,1]{\rm cos}\,\iota\in[-1,1], ψ∈[0,360∘]\psi\in[0,360^{\circ}], and ψc∈[0,360∘]\psi_{\rm c}\in[0,360^{\circ}]. For the NS’s mass, we adopt the mass distribution based on the analysis of the latest O3 observation (figure 7 of ref. [111]). Without loss of generality, the merger time is chosen to tc=0t_{\rm c}=0 in our analysis. The simulated GW dark sirens satisfy the redshift distributions corresponding to ET, CE, 2CE, and ET2CE shown in figure 2.

When using BNS tidal deformation measurements to obtain redshift information, the GW tidal phases need to be precisely measured to obtain the source-frame mass information (related to BB and CC). Combining with the observed chirp mass information obtained from the measurements of observed GW frequency ff and its time derivative f˙\dot{f}, the redshift information could be obtained by breaking the parameter degeneracies between mass and redshift. In our analysis, we consider the correlation between the twelve parameters to calculate the measurement errors of redshifts. Using eq. (2.12), we can obtain the measurement errors of redshifts, Δ​z=(F−1)22\Delta z=\sqrt{(F^{-1})_{22}}.

In the upper and middle panels of figure 3, we show the distributions of Δ​z/z\Delta z/z for ET, CE, 2CE, and ET2CE. In order to be more intuitive, in the lower panel of figure 3, we use the interpolation method to fit the histograms above. We see that the measurement precisions of redshifts are mainly 5%5\%–20%20\%. ET2CE gives the best redshift measurement, followed by 2CE, CE, and ET. The ET case shows evident difference compared with the other three cases, while the CE, 2CE, and ET2CE cases also have some slight differences.

In figure 4, we show the simulated GW standard sirens observed by ET, CE, 2CE, and ET2CE. For simplicity, we show only 30 data points for each observation. To account for fluctuations in the measured values resulting from actual observations, we represent the standard siren data points with Gaussian randomization. Specifically, the central value is populated according to a Gaussian distribution with mean equal to the fiducial value and standard deviation equal to the measurement error.

Refer to caption
Figure 4: The simulated standard sirens observed by ET, CE, 2CE, and ET2CE. Here we show only 30 data points for each observation for simplicity. We actually simulate 545763, 1017618, 1317183, and 1815501 standard sirens for ET, CE, 2CE, and ET2CE.

2.4 Method of constraining cosmological parameters

Given a set of NN GW data 𝒟→\vec{\mathcal{D}}, the posterior probability distribution of the cosmological parameters Ω→\vec{\Omega} is

p​(Ω→|𝒟→)∝ℒ​(𝒟→|Ω→)​p​(Ω→),p(\vec{\Omega}|\vec{\mathcal{D}})\propto\mathcal{L}(\vec{\mathcal{D}}|\vec{\Omega})p(\vec{\Omega}), (2.15)

where p​(Ω→)p(\vec{\Omega}) is the prior on cosmological parameters Ω→\vec{\Omega}. The likelihood function for the GW data ℒ​(𝒟→|Ω→)\mathcal{L}(\vec{\mathcal{D}}|\vec{\Omega}) is given by

ℒ​(𝒟→|Ω→)=∏iNp​(dL,obsi|dL​(zi,Ω→))​p​(zobsi|zi)​p​(zi),\mathcal{L}(\vec{\mathcal{D}}|\vec{\Omega})=\prod_{i}^{N}p(d_{\mathrm{L,obs}}^{i}|d_{\mathrm{L}}(z^{i},\vec{\Omega}))p(z_{\mathrm{obs}}^{i}|z^{i})p(z^{i}), (2.16)

where ii represents the ii-th GW event, p​(dL,obsi|dL​(zi,Ω→))p(d_{\mathrm{L,obs}}^{i}|d_{\mathrm{L}}(z^{i},\vec{\Omega})) and p​(zobsi|zi)p(z_{\mathrm{obs}}^{i}|z^{i}) are two Gaussian distributions with the standard deviations calculated from eq. (2.12), and p​(zi)p(z^{i}) is the redshift prior that is assumed to be uniform for simplicity. For the dark siren method, the redshifts zz of the GW events are not precisely measured. Therefore, the redshift zz should be marginalized out to obtain the constraints on cosmological parameters, i.e., eq. (2.16) should be rewritten as

ℒ​(𝒟→|Ω→)=∏iN∫p​(dL,obsi|dL​(zi,Ω→))​p​(zobsi|zi)​dzi.\mathcal{L}(\vec{\mathcal{D}}|\vec{\Omega})=\prod_{i}^{N}\int p(d_{\mathrm{L,obs}}^{i}|d_{\mathrm{L}}(z^{i},\vec{\Omega}))p(z_{\mathrm{obs}}^{i}|z^{i})\mathrm{d}z^{i}. (2.17)

Note that in the above equation p​(zi)p(z^{i}) is not explicitly shown because it is assumed to be a constant.

Here, we also note that the method described above is valid under two assumptions: (i) the measurements of zz and dLd_{\rm L} are assumed to be nearly independent, with dLd_{\rm L} measured from the GW amplitude and zz measured from the GW phase, and (ii) the detection probability (selection bias) does not depend on the observed values of zz and dLd_{\rm L}. In a more realistic approach, the observed values of zz and dLd_{\rm L} may deviate from the simulated values. In such cases, it is necessary to use the observed data to calculate SNRs and the detection probability. If assumption (ii) is not satisfied, the Bayesian inference could be corrected by adding the redshift prior conditioned on detection, which is basically the detection rate as a function of redshift (see ref. [131] for an example).

Table 2: The absolute errors (1σ\sigma) and the relative errors of the cosmological parameters in the wwCDM and w0​waw_{0}w_{a}CDM models by using the mock data of ET, CE, 2CE, and ET2CE.

Model Parameter ET CE 2CE ET2CE wwCDM σ​(Ωm)\sigma(\Omega_{m}) 0.00550.0055 0.00280.0028 0.00130.0013 0.00110.0011 σ​(H0)\sigma(H_{0}) 0.2700.270 0.1870.187 0.1250.125 0.1030.103 σ​(w)\sigma(w) 0.0340.034 0.0200.020 0.0120.012 0.0100.010 ε​(Ωm)\varepsilon(\Omega_{\rm m}) 1.75%1.75\% 0.89%0.89\% 0.43%0.43\% 0.32%0.32\% ε​(H0)\varepsilon(H_{0}) 0.40%0.40\% 0.27%0.27\% 0.18%0.18\% 0.15%0.15\% ε​(w)\varepsilon(w) 3.39%3.39\% 2.01%2.01\% 1.19%1.19\% 0.95%0.95\% w0​waw_{0}w_{a}CDM σ​(Ωm)\sigma(\Omega_{m}) 0.02550.0255 0.01130.0113 0.00590.0059 0.00460.0046 σ​(H0)\sigma(H_{0}) 0.3660.366 0.2500.250 0.1850.185 0.1600.160 σ​(w0)\sigma(w_{0}) 0.0500.050 0.0340.034 0.0240.024 0.0200.020 σ​(wa)\sigma(w_{a}) 0.5040.504 0.2710.271 0.1570.157 0.1260.126 ε​(Ωm)\varepsilon(\Omega_{\rm m}) 8.09%8.09\% 3.56%3.56\% 1.89%1.89\% 1.48%1.48\% ε​(H0)\varepsilon(H_{0}) 0.54%0.54\% 0.37%0.37\% 0.27%0.27\% 0.24%0.24\% ε​(w0)\varepsilon(w_{0}) 5.01%5.01\% 3.42%3.42\% 2.39%2.39\% 2.04%2.04\%

3 Results

Refer to caption
Refer to caption
Figure 5: Two-dimensional marginalized contours (68.3%68.3\% and 95.4%95.4\% confidence level) in the Ωm\Omega_{\rm m}–H0H_{0} and ww–Ωm\Omega_{\rm m} planes using the mock GW dark siren data of ET, CE, 2CE, and ET2CE. Here, the dotted lines indicate the fiducial values of cosmological parameters preset in the simulation.

We use the simulated GW dark siren data from ET, CE, 2CE, and ET2CE and by performing the MCMC analysis [132] we use the mock data of dark sirens to constrain the wwCDM and w0​waw_{0}w_{a}CDM models. The constraint results are shown in figures 5 and 6 and summarized in table 2. For a parameter ξ\xi, we use σ​(ξ)\sigma(\xi) and ε​(ξ)\varepsilon(\xi) to represent its absolute (1σ\sigma) and relative errors, respectively, with ε​(ξ)\varepsilon(\xi) defined as σ​(ξ)/ξ\sigma(\xi)/\xi.

We first take a look at the wwCDM model. In the left panel of figure 5, we show the constraints in the Ωm\Omega_{\rm m}–H0H_{0} plane. We can clearly see that GWs place quite tight constraints on the parameters. Owing to the fact that the detection sensitivity of CE is better than that of ET, ET shows the worst constraint results among the four cases. However, even so, ET gives σ​(Ωm)=0.0055\sigma(\Omega_{\rm m})=0.0055, σ​(H0)=0.270\sigma(H_{0})=0.270, and σ​(w)=0.034\sigma(w)=0.034 which are comparable with the results of σ​(Ωm)=0.0065\sigma(\Omega_{\rm m})=0.0065, σ​(H0)=0.660\sigma(H_{0})=0.660, and σ​(w)=0.028\sigma(w)=0.028 by the latest Planck TT,TE,EE+lowE+BAO+Pantheon data [32]. Meanwhile, ET and CE show positive correlation between Ωm\Omega_{\rm m} and H0H_{0}, while there is almost no correlation between Ωm\Omega_{\rm m} and H0H_{0} when using the mock GW dark siren data of 2CE and ET2CE. However, the obvious anti-correlation is obtained when using the traditional EM cosmological probes. In the right panel of figure 5, we show the constraints in the ww–Ωm\Omega_{\rm m} plane. We find that the 2CE network can give quite tight constraints on the cosmological parameters with the constraint precisions of Ωm\Omega_{\rm m}, H0H_{0}, and ww being 0.43%0.43\%, 0.18%0.18\%, and 1.19%1.19\%, respectively. So, even for the EoS of dark energy, the standard of precision cosmology is nearly achieved. Furthermore, we find that using the ET2CE mock data, the constraint precisions of Ωm\Omega_{\rm m}, H0H_{0}, and ww are 0.32%0.32\%, 0.15%0.15\%, and 0.95%0.95\%, respectively, all meeting the standard of precision cosmology.

In figure 6, we show the constraint results of the w0​waw_{0}w_{a}CDM model. Among the four cases, ET gives the worst constraint results. CE gives much better constraints than ET. 2CE gives similar constraint results with ET2CE. Quantitatively, ET gives σ​(w0)=0.050\sigma(w_{0})=0.050 and σ​(wa)=0.504\sigma(w_{a})=0.504, comparable with the results of σ​(w0)=0.064\sigma(w_{0})=0.064 and σ​(wa)=0.300\sigma(w_{a})=0.300 by the latest Planck TT,TE,EE+lowE+BAO+Pantheon data [32]. Meanwhile, the constraint results of ET are also basically consistent with those given in ref. [95]. For the case of using CE, the constraint results are much better than those of ET, also better than those of the latest Planck TT,TE,EE+lowE+BAO+Pantheon data. What’s more, when considering 2CE and ET2CE, the constraint precisions of w0w_{0} are 2.39%2.39\% (2CE) and 2.04%2.04\% (ET2CE). For waw_{a}, ET2CE gives σ​(wa)=0.126\sigma(w_{a})=0.126.

Refer to caption
Refer to caption
Figure 6: Two-dimensional marginalized contours (68.3%68.3\% and 95.4%95.4\% confidence level) in the Ωm\Omega_{\rm m}–H0H_{0} and w0w_{0}–waw_{a} planes using the mock GW dark siren data of ET, CE, 2CE, and ET2CE. Here, the dotted lines indicate the fiducial values of cosmological parameters preset in the simulation.

4 Discussion

In section 3, we compare the constraint results of the GW dark sirens (tidal measurements) with the current cosmological observations, from which the potential of the GW dark sirens can be clearly seen. Since third-generation GW detectors will begin observation in the 2030s, it is also necessary to show the comparisons with other cosmological observations in the same period.

4.1 Comparison with future other GW standard siren observations

In the 3G GW detector network era, lots of bright sirens could also be used to study the cosmic expansion history. Belgacem et al. [90] used the mock bright siren data (based on the 10-year observation of ET2CE and THESEUS) to measure dark energy and obtained σ​(w)=0.034\sigma(w)=0.034 for the optimistic case of observing EM counterparts. When combining the mock bright siren data of ET2CE with the current CMB+BAO+SN data, the constraint errors could be reduced to σ​(w)=0.013\sigma(w)=0.013 that is comparable with the constraint results of 2CE in this work. For the case of the w0​waw_{0}w_{a}CDM model, the combination of the mock bright siren data and the CMB+BAO+SN data gives σ​(w0)=0.025\sigma(w_{0})=0.025 and σ​(wa)=0.137\sigma(w_{a})=0.137 [90] that are also comparable with those of 2CE in this work.

In the same period of the 3G GW detectors, the space-based GW detectors will begin detecting GWs generated by the massive or supermassive black hole binaries (SMBHBs). Wang et al. [50] used the mock bright siren data of the LISA-Taiji network to measure dark energy. Zhu et al. [133] used the mock bright siren data of LISA-TianQin network to measure dark energy. However, owing to the fact that the number of the detectable bright sirens of SMBHBs are small (about 𝒪​(10)\mathcal{O}(10)) within 5-year observation, the constraint results are all weaker than those of this work.

4.2 Comparison with future EM surveys

In the next decades, the fourth-generation dark energy surveys will also precisely measure dark energy [134, 135]. Pierel et al. [134] forecasted the ability of cosmological parameter estimation using strongly lensed SN from the Roman Space Telescope (previously WFIRST). The CMB+BAO+SN+UltraDeep lensed SN data give σ​(w)=0.01\sigma(w)=0.01 in the wwCDM model, and σ​(w0)=0.02\sigma(w_{0})=0.02 and σ​(wa)=0.14\sigma(w_{a})=0.14 in the w0​waw_{0}w_{a}CDM model, comparable with those of 2CE in this work. Recently, Euclid Collaboration [135] forecasted the constraint results of Euclid. The pessimistic case of Euclid (see ref. [135] for more details) gives σ​(w0)=0.038\sigma(w_{0})=0.038 and σ​(wa)=0.14\sigma(w_{a})=0.14 that are comparable with those of 2CE. If considering the combination of CMB-S4 and Euclid (pessimistic), the constraint results are worse than those of ET2CE. While considering the optimistic case of Euclid (see ref. [135] for more details), the absolute errors are reduced to σ​(w0)=0.021\sigma(w_{0})=0.021 and σ​(wa)=0.073\sigma(w_{a})=0.073 that are comparable with those of ET2CE. If considering the combination of CMB-S4 and Euclid (optimistic), the constraint results are also comparable with those of ET2CE.

4.3 Comparison with future other cosmological probes

In the future, neutral hydrogen 21 cm intensity mapping (IM) will become an important tool to explore the nature of dark energy [136, 137, 138, 139, 140, 141, 142]. Bull et al. [136] obtained σ​(w0)=0.035\sigma(w_{0})=0.035 and σ​(wa)=0.163\sigma(w_{a})=0.163 by using the SKA1-MID+MeerKAT data, which are comparable with those of 2CE. Xu et al. [137] obtained σ​(w0)=0.082\sigma(w_{0})=0.082 and σ​(wa)=0.21\sigma(w_{a})=0.21 by using the full-scale Tianlai data, which are comparable with those of ET. Wu et al. [56] forecasted cosmological parameters using the combination of four promising cosmological probes, 21 cm IM, GW, fast radio burst, and strong gravitational lensing, and obtained σ​(w)=0.020\sigma(w)=0.020, comparable with that of CE.

5 Conclusion

In the next decades, the third-generation GW detectors will detect large amounts of (𝒪​(105)\mathcal{O}(10^{5}) per year) GW events generated by the BNS mergers. By measuring the additional tidal deformation-phase term, the redshifts of the GW sources could be obtained. Hence, the cosmic expansion history can be investigated using solely GW observations. In this work, based on the three-year observation, we simulate such GW dark sirens of ET, CE, 2CE, and ET2CE to perform cosmological analysis. Note that the assumption of the EoS of NS being perfectly known is adopted in this work. We wish to investigate whether the Hubble constant and the EoS of dark energy can both be precisely measured using only GW dark sirens detected by the 3G GW observations.

We find that ET gives the worst constraints among the four cases. Even so, ET gives σ​(w)=0.034\sigma(w)=0.034 that is comparable with the latest constraint result of Planck TT,TE,EE+lowE+ BAO + Pantheon. Meanwhile, the constraint precision of ww using the mock 2CE data is close to 1%1\%, which almost meets the standard of precision cosmology. The 3G GW detector network has a strong capability to constrain the EoS parameter of dark energy ww, with a precision of 0.95%. For the Hubble constant, the single ET observatory gives quite tight constraint with the precision being 0.40%0.40\%. Using the ET2CE network, we obtain the constraint precision of H0H_{0} being 0.15%0.15\%.

We also discuss the case of the two-parameter dark-energy model. We find that ET gives comparable constraint results with those of the latest Planck TT,TE,EE+lowE+ BAO + Pantheon data. CE gives much better constraint results than those of ET. ET2CE gives slightly better constraint results than 2CE. Using ET2CE, we obtain the constraint precision of w0w_{0} being 2.04%2.04\%. Moreover, ET2CE gives σ​(wa)=0.126\sigma(w_{a})=0.126 that is much better than the result of the latest Planck TT,TE,EE+lowE+BAO+Pantheon data.

In addition, we discuss the comparison of the constraint results of GW dark sirens (tidal measurements) with those of the same-period other cosmological observations. We find that GW dark sirens show strong ability to measure the Hubble constant and dark energy. The constraints on the Hubble constant and the EoS of dark energy using GW dark sirens of ET2CE are comparable with those of the combination of the fourth-generation dark energy surveys (optimistic case of Euclid) and CMB-S4.

Our results show that the precision cosmology can be achieved by solely using the GW dark sirens observed by the 3G GW detectors. It can be expected that the 3G GW detectors would play a key role in helping solve the Hubble tension and reveal the fundamental nature of dark energy.

Acknowledgments

We thank Ji-Yu Song, Jiming Yu, Liang-Gui Zhu, Arnab Dhani, Ssohrab Borhanian, Peng-Ju Wu, Ling-Feng Wang, Tao Han, Jing-Zhao Qi, Bo Wang, and Rui-Qi Zhu for helpful discussions. This work was supported by the National SKA Program of China (Grants Nos. 2022SKA0110200 and 2022SKA0110203), the National Natural Science Foundation of China (Grants Nos. 11975072, 11875102, and 11835009), and the 111 Project (Grant No. B16009).

References