Multiple Measurements of Gravitational Waves Acting as Standard Probes: Model-independent Constraints on the Cosmic Curvature with DECIGO

The spatial curvature parameter Ω_k is an important quantity that is related to many fundamental issues in modern cosmology, such as the structure and evolution of the universe (Ichikawa & Takahashi 2006; Zolnierowski & Blanchard 2015; Cao et al. 2019a; Qi et al. 2019a). Specifically, the study of the cosmic curvature can effectively test the fundamental assumption of modern cosmology that the universe is homogeneous and isotropic and is well described by the FLRW metric. The most popular theory of the very early universe proposes that our universe once went through an exponential phase of expansion, which indicates that the radius of curvature of the universe should be very large and the cosmic curvature should be close to zero (Weinberg et al. 2013). Current cosmological observations, particularly the latest Planck2018 results, which combined the cosmic microwave background and baryon acoustic oscillation (BAO) measurements, strongly favor a flat universe: Ω_k = 0.0007 ± 0.0019 (Planck Collaboration et al. 2020) (TT,TE,EE+lowE+lensing+BAO). However, the combination of the Planck2018 TT,TE,EE+lowE power spectra data alone marginally favors a mildly closed universe: ${{\rm{\Omega }}}_{k}=-{0.044}_{-0.015}^{+0.018}$ (Planck Collaboration et al. 2020; Di Valentino et al. 2020). In addition, such a stringent constraint on the cosmic curvature relies heavily on the assumption of a specific cosmological model (the cosmological constant plus cold dark matter model, i.e., the ΛCDM model). However, recent analysis indicated that the flat-universe assumption may lead to an incorrect reconstruction of the dark energy equation of state and cause tension between the ΛCDM and dynamical dark energy models (Ichikawa et al. 2006; Clarkson et al. 2007; Gong & Wang 2007; Virey et al. 2008; Li et al. 2018; Cao et al. 2019b). The strong degeneracy between the cosmic curvature Ω_k and the dark energy equation of state w makes it difficult to constrain these two parameters simultaneously in a non-flat wCDM model. Therefore, it would be better to measure spatial curvature in geometric and model-independent ways.

The distance sum rule has been proposed as a model-independent method to constrain the curvature of the universe, and is generally implemented using strong lensing observations (Cao & Zhu 2012a; Cao et al. 2012b, 2015; Ma et al. 2019a) with other distance measurements, such as Type Ia supernovae (SNe Ia) (Räsänen et al. 2015; Denissenya et al. 2018) or quasars (QSOs) (Qi et al. 2019b; Zhou & Li 2020; Lian et al. 2021). Another model-independent way to constrain the cosmic curvature is by comparing the theoretical comoving distance, which is inferred from the observational Hubble parameter data (OHD), with the observed luminosity distances D_L (z) (or the angular diameter distance D_A (z)) (Clarkson et al. 2008; Shafieloo & Clarkson 2010; Mortsell & Onsson 2011; Sapone et al. 2014; Cai et al. 2016). This test has been fully implemented with various observational data. For example, Wang et al. (2020) recently used the latest Pantheon SNe Ia data combined with cosmic chronometers (CCs) to constrain the cosmic curvature; a well-measured quasar sample (Risaliti & Lusso 2015, 2017; Lusso & Risaliti 2016, 2017) also shows great potential for probing the cosmic curvature (Wei & Melia 2020); Cao et al. (2019b) proposed an improved model-independent test of cosmic curvature with ultracompact structures in radio quasars as standard rulers, which shows no evidence for a deviation from a flat universe; and Takada & Doré (2015) proposed that the combination of radial and angular diameter distances from future BAO experiments can be used for studying the curvature parameter. One can also use an analytical equation to reconstruct Ω_k at different redshifts (Clarkson et al. 2007), which also requires observational data on the Hubble parameter H(z) and luminosity distance D_L (z), plus the first derivative of the latter. Such a method helps us study the evolution of the spatial curvature and provides a direct geometric way to test the assumption of cosmic homogeneity and isotropy. Ω_k (z) reconstructed in this way, using current observational data, the combination of cosmic chronometer data with 1598 quasars (Risaliti & Lusso 2018), and the Pantheon catalog of 1048 SNe Ia (Scolnic et al. 2018), shows good agreement with a flat universe at different redshifts (Liu et al. 2020).

On the other hand, gravitational-wave (GW) observations soon caught people's attention with the discovery of the first GW event, GW150914 (Abbott et al. 2016). As standard sirens, GW signals from inspiraling and merging compact binaries encode distance information (Schutz 1986) and provide an absolute measurement of the luminosity distance. Only if these binary mergers are accompanied by short–hard γ-ray bursts (shGRBs) can they be observed through both electromagnetic (EM) and gravitational waves. The joint detection of GW170817 (Abbott et al. 2017) has detected its electromagnetic counterparts from the merger of binary neutron stars (NSs). Knowing the redshifts of the sources, these GW signals can be used for cosmology. Also, the accelerating expansion of the universe would cause an additional phase shift in the gravitational waveform, which allows us to measure the cosmic acceleration (or redshift drift) directly (Cutler & Holz 2009; Nishizawa et al. 2011). Thus, it is possible for us to break the degeneracy and measure cosmic curvature independently in the GW domain.

However, the measurement of the cosmic acceleration requires a high-precision GW detection, particularly at lower frequencies when the binary remains in its inspiraling phase. For GW signals from a neutron star binary at a redshift of 1, the acceleration of the universe's expansion would cause a phase delay of only 1 s during the 10 yr probe (Seto et al. 2001; Nishizawa et al. 2012). As a future spaceborne GW detector, the Decihertz Interferometer Gravitational-wave Observatory (DECIGO) (Seto et al. 2001; Kawamura et al. 2011) is designed to improve the detection sensitivity of GWs at lower frequencies, with its most sensitive frequencies between 0.1 and 10 Hz. DECIGO will have four clusters of spacecraft, each cluster consisting of three spacecraft with three Fabry–Perot Michelson interferometers, whose arm length is 1000 km, improving the determination of their position in the sky. The expected sensitivity of DECIGO is 10⁻²⁵ Hz^−1/2 for two clusters at the same position for three years of mission, which enables the early detection of inspiraling sources. Therefore, DECIGO would create an unprecedented opportunity to precisely measure the cosmic acceleration from GW signals and make GWs a more precise standard siren.

Due to the lack of GW events with observed EM counterparts, we simulated 10,000 GW events from DECIGO within the redshift range 0–5. We applied two model-independent methods to measure the cosmic curvature: numerical constraint and reconstruction. For comparison, we also use Pantheon SNe Ia and OHD to estimate Ω_k . The remainder of this paper is organized as follows. In Section 2, we introduce the simulated GW data and methodology used in this study. In Section 3, we present the results of our study and provide some discussion by comparison. Finally, the general conclusions are summarized in Section 4.

Assuming that the universe is homogeneous and isotropic on large scales, it can be described by the FLRW metric:

$$\begin{eqnarray}&&{{ds}}^{2}=-{{dt}}^{2}+{a}^{2}(t)\left[\displaystyle \frac{{{dr}}^{2}}{1-{{Kr}}^{2}}+{r}^{2}(d{\theta }^{2}+{\sin }^{2}\theta d{\phi }^{2})\right],\end{eqnarray} \tag{ 1 }$$

where t is the cosmic time and (r, θ, ϕ) are the comoving spatial coordinates. The scale factor a(t) is the only gravitational degree of freedom and its evolution is determined by the matter and energy of the universe. The dimensionless curvature K = −1, 0, or +1 corresponds to open, flat, and closed universes respectively. With such a metric, the luminosity distance D_L (z) can be expressed as

$$\begin{eqnarray}{D}_{L}(z)=\left\{\begin{array}{lll}\displaystyle \frac{c(1+z)}{{H}_{0}\sqrt{| {{\rm{\Omega }}}_{{\rm{k}}}| }}\sinh \left[\sqrt{| {{\rm{\Omega }}}_{{\rm{k}}}| }{\displaystyle \int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )}\right] & \mathrm{for} & {{\rm{\Omega }}}_{k}\gt 0,\\ \displaystyle \frac{c(1+z)}{{H}_{0}}{\displaystyle \int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )} & \mathrm{for} & {{\rm{\Omega }}}_{K}=0,\\ \displaystyle \frac{c(1+z)}{{H}_{0}\sqrt{| {{\rm{\Omega }}}_{{\rm{k}}}| }}\sin \left[\sqrt{| {{\rm{\Omega }}}_{{\rm{k}}}| }{\displaystyle \int }_{0}^{z}\displaystyle \frac{{dz}^{\prime} }{E(z^{\prime} )}\right] & \mathrm{for} & {{\rm{\Omega }}}_{k}\lt 0.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

The dimensionless Hubble parameter E(z) is defined as H(z)/H₀, where H(z) is the expansion rate of the universe and H₀ is the Hubble constant. The curvature parameter Ω_k is related to K as ${{\rm{\Omega }}}_{k}=-{c}^{2}K/{\left({a}_{0}{H}_{0}\right)}^{2}$, where c is the speed of light.

2.1. Measuring the Distance and the Acceleration Parameter with GWs

As a standard siren, a GW signal allows us to measure the luminosity distance directly without relying on the cosmic distance ladder. Concurrently, if the expansion of the universe is accelerating, we might find an additional phase shift in the gravitational waveform, from which we can measure the acceleration parameter X, and thus obtain the Hubble parameter. Therefore, both the luminosity distance and the Hubble parameter can be obtained from a single GW signal. To apply this D_L –X relation to cosmological studies, we still need information about the redshifts, which requires the GW sources to be neutron star binaries (NS–NS) or black hole–neutron star binaries (BH–NS), the origins of kilonovae or shGRBs. As long as the corresponding EM counterparts are observed from GW events, we can obtain their redshifts and apply them to cosmology (discussed in Section 2.2).

In this study, we focus on the GW signals from NS–NS binary systems with component masses m₁ and m₂. One can define the chirp mass M_c = M η^3/5, and the redshifted chirp mass ${{ \mathcal M }}_{z}\equiv M(1+{z}_{c}){\eta }^{3/5}$, where M = m₁ + m₂ is the total mass of the binary system, η = m₁ m₂/M² is the symmetric mass ratio, and z_c is the source redshift at coalescence.

We first derive the correction to the GW phase due to the accelerating expansion of the universe. The observed gravitational waveform can be represented by h(Δt), where Δt ≡ t_c − t denotes the time to coalescence measured in the observer frame, with t_c representing the coalescence time. The Fourier component of this waveform can be written as

$$\begin{eqnarray}&&\tilde{h}(f)={\int }_{-\infty }^{\infty }{e}^{2\pi {ift}}h({\rm{\Delta }}t){dt}.\end{eqnarray} \tag{ 3 }$$

Then, we define Δt_e as the time to coalescence measured in the source frame, and ΔT ≡ (1 + z_c )Δt_e as the redshifted coalescence time. The coalescence times in these two different frames have a relation: ⁴

$$\begin{eqnarray}&&{\rm{\Delta }}t={\rm{\Delta }}T+X\left({z}_{c}\right){\rm{\Delta }}{T}^{2},\end{eqnarray} \tag{ 4 }$$

which is measured by the acceleration parameter X(z) (Seto et al. 2001; Takahashi & Nakamura 2005) and the correction term due to cosmic acceleration. The acceleration parameter X(z) is defined as ⁵

$$\begin{eqnarray}&&X(z)\equiv \displaystyle \frac{{H}_{0}}{2}\left(1-\displaystyle \frac{H(z)}{(1+z){H}_{0}}\right).\end{eqnarray} \tag{ 5 }$$

By substituting Equation (4) into Equation (3) and applying the stationary phase approximation (Cutler & Flanagan 1994), we can transform the Fourier component of the waveform into the frequency domain:

$$\begin{eqnarray}&&\tilde{h}(f)={\left.{e}^{i{{\rm{\Psi }}}_{\mathrm{acc}}(f)}\tilde{h}(f)\right|}_{\mathrm{no}\ \mathrm{acc}},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Psi }}}_{\mathrm{acc}}(f) & \equiv & -2\pi {fX}\left({z}_{c}\right){\rm{\Delta }}T{\left(f\right)}^{2}\\ & = & -{{\rm{\Psi }}}_{N}(f)\displaystyle \frac{25}{768}X\left({z}_{c}\right){{ \mathcal M }}_{z}{x}^{-4}\end{array}\end{eqnarray} \tag{ 7 }$$

corresponds to the gravitational waveform with cosmic acceleration, with $x\equiv {\left(\pi {{ \mathcal M }}_{z}f\right)}^{2/3}$ and ${{\rm{\Psi }}}_{N}(f)\,\equiv \tfrac{3}{128}{\left(\pi {{ \mathcal M }}_{z}f\right)}^{-5/3}$. The gravitational waveform without cosmic acceleration is given by

$$\begin{eqnarray}&&\tilde{h}(f)\left|{}_{\mathrm{no}\,\mathrm{acc}}\right.=\displaystyle \frac{\sqrt{3}}{2}{ \mathcal A }{f}^{-7/6}{e}^{i{\rm{\Psi }}(f)}\left[\displaystyle \frac{5}{4}{A}_{\mathrm{pol},\alpha }(t(f))\right]{e}^{-i\left({\varphi }_{\mathrm{pol},\alpha }+{\varphi }_{D}\right)},\end{eqnarray} \tag{ 8 }$$

with the amplitude written as

$$\begin{eqnarray}&&{ \mathcal A }=\displaystyle \frac{1}{\sqrt{30}{\pi }^{2/3}}\displaystyle \frac{{{{ \mathcal M }}_{z}}^{5/6}}{{D}_{L}}.\end{eqnarray} \tag{ 9 }$$

The polarization amplitude A_pol,α (t), the polarization phases φ_pol,α (t) (α = I, II represents the number of individual detectors ⁶ ), and the Doppler phase φ_D (t) are given in Yagi & Tanaka (2010). For the phase of Ψ(f), we use the restricted-2PN (PN = post-Newtonian) waveform including spin–orbit coupling at the 1.5PN order ⁷ (Kidder et al. 1993), which is also reported by Yagi & Tanaka (2010). There, we can extract the acceleration parameter X and luminosity distance D_L from the phase and amplitude of the GWs, respectively.

According to the waveform in Equation (6), we take the binary parameters as

$$\begin{eqnarray}&&{\theta }^{i}=\left(\mathrm{ln}{{ \mathcal M }}_{z},\mathrm{ln}\eta ,\beta ,{t}_{c},{\phi }_{c},{\bar{\theta }}_{{\rm{S}}},{\bar{\phi }}_{{\rm{S}}},{\bar{\theta }}_{{\rm{L}}},{\bar{\phi }}_{{\rm{L}}},{D}_{L},X\right).\end{eqnarray} \tag{ 10 }$$

where β is the spin–orbit coupling parameter and ϕ_c represents the coalescence phase. $({\bar{\theta }}_{{\rm{S}}},{\bar{\phi }}_{{\rm{S}}})$ indicates the direction of the source in the solar barycentric frame and $({\bar{\theta }}_{{\rm{L}}},{\bar{\phi }}_{{\rm{L}}})$ specifies the direction of the orbital angular momentum.

One can use Fisher analysis to estimate the measurement accuracies of the binary parameters θⁱ . The measurement accuracy is given by ${\rm{\Delta }}{\theta }^{i}\equiv \sqrt{{\left({{\rm{\Gamma }}}^{-1}\right)}_{{ii}}}$ (Cutler & Flanagan 1994), with the Fisher matrix:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}=4{\rm{Re}}{\int }_{{f}_{min}}^{{f}_{max}}\displaystyle \frac{{{\rm{\partial }}}_{i}\tilde{h}\ast (f){{\rm{\partial }}}_{j}\tilde{h}(f)}{{S}_{h}(f)}{df},\end{eqnarray} \tag{ 11 }$$

where S_h (f) is the analytical expression of the DECIGO noise power spectrum (Kawamura et al. 2006; Yagi & Seto 2011; Kawamura et al. 2019). The lower cutoff of frequency ${f}_{min}\,={\left(256/5\right)}^{-3/8}{\pi }^{-1}{{ \mathcal M }}_{z}^{-5/8}{\rm{\Delta }}{T}_{{\rm{obs}}}^{-3/8}$ corresponds to the frequency at which coalescence begins to be observed, with ΔT_obs representing the observation time. f_max is the higher cutoff frequency of the detector and is set equal to 100 Hz. In addition, for DECIGO there will be eight uncorrelated interferometric signals (Kawamura et al. 2011; Yagi & Seto 2011); thus, the Fisher matrix above should be multiplied by a factor of 8.

For the convenience of the above calculation, we set m₁ = m₂ = 1.4 M_⊙ and take t_c = ϕ_c = β = 0. For each fiducial redshift z_c , we randomly generate 10⁴ sets of $({\bar{\theta }}_{{\rm{S}}},{\bar{\phi }}_{{\rm{S}}},{\bar{\theta }}_{{\rm{L}}},{\bar{\phi }}_{{\rm{L}}})$, and for each set, we calculate the uncertainty $\sqrt{{\left({{\rm{\Gamma }}}^{-1}\right)}_{{ii}}}$. Therefore, only two parameters need to be estimated (D_L and X). By marginalizing other parameters in the Fisher matrix, this calculation converts into a two-dimensional submatrix of D_L and X, with the instrumental uncertainty of X being

$$\begin{eqnarray}&&{\sigma }_{X}={8}^{-1/2}\left[{\left({{\rm{\Gamma }}}^{-1}\right)}_{{jj}}^{1/2}\right].\end{eqnarray} \tag{ 12 }$$

Similarly, the measurement error of luminosity distance is estimated as

$$\begin{eqnarray}&&{\sigma }_{{D}_{L}}^{\mathrm{instr}}={8}^{-1/2}\left[{\left({{\rm{\Gamma }}}^{-1}\right)}_{{ii}}^{1/2}\right].\end{eqnarray} \tag{ 13 }$$

Considering the lensing uncertainty caused by the weak lensing effect ${\sigma }_{{D}_{L}}^{\mathrm{lens}}=0.05{{zD}}_{L}$ (Sathyaprakash et al. 2010), the luminosity distance error is considered to be

$$\begin{eqnarray}&&{\sigma }_{{D}_{L}}=\sqrt{{\left({\sigma }_{{D}_{L}}^{\mathrm{instr}}\right)}^{2}+{\left({\sigma }_{{D}_{L}}^{\mathrm{lens}}\right)}^{2}}.\end{eqnarray} \tag{ 14 }$$

We sample the redshifts of GW sources from the merger rate of double compact objects, which reflects the star formation history (Dominik et al. 2013), taking the data from the so-called "rest frame rates" in the cosmological scenario. ⁸ In our simulation, we assume a flat ΛCDM model as the proposed fiducial model with the cosmological parameters derived from Planck2018 measurements (Planck Collaboration et al. 2020). The central values of the simulated 10,000 X(z) and D_L (z) measurements are shown in Figure 1, and the acceleration parameters and their corresponding errors are shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.2. Redshift Determination from Optical Follow-up Observations

2.3. Numerical Constraints on Cosmic Curvature Ω_k (z)

GWs can provide direct measurements of the luminosity distance D_L (z) and the cosmic acceleration X(z), as we simulated in Section 2.1. Then, one can simply confront the distance D_L (z) with theoretical distances ${D}_{L}^{X}(z)$ inferred from Equation (2), where curvature parameter Ω_k is considered and E(z) can be derived from the definition equation of X(z) in Equation (5), to provide a numerical constraint on the cosmic curvature Ω_k (z). This process could be achieved by minimizing the χ² statistic

$$\begin{eqnarray}&&{\chi }^{2}({{\rm{\Omega }}}_{k})=\sum _{i=1}^{10,000}\displaystyle \frac{{\left[{D}_{L,i}^{X}({{\rm{\Omega }}}_{k})-{D}_{L,i}\right]}^{2}}{{\sigma }_{{D}_{L}^{X},i}^{2}+{\sigma }_{{D}_{L},i}^{2}},\end{eqnarray} \tag{ 15 }$$

where ${\sigma }_{{D}_{L}^{X}}$ can be derived from Equation (2) and ${\sigma }_{{D}_{L}}$ is given in Equation (14). We use the Markov Chain Monte Carlo method (Foreman-Mackey et al. 2013) to obtain the best-fit value of Ω_k and its uncertainty, which is a model-independent estimation of the cosmic curvature.

2.4. Individual Measurements of Cosmic Curvature Ω_k (z)

Future space-based GW detectors, such as DECIGO, will aim to detect a large number of NS binaries at much higher redshifts, which will enable the reconstruction of cosmic curvature in the early universe. The cosmic curvature Ω_k can be expressed as (Clarkson et al. 2007)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{k}=\displaystyle \frac{{\left[H(z){D}^{{\prime} }(z)\right]}^{2}-{c}^{2}}{{\left[{H}_{0}D(z)\right]}^{2}},\end{eqnarray} \tag{ 16 }$$

where the expansion rate of the universe H(z) can be obtained from the GW measurement of X(z), the transverse comoving distance D(z) is simply related to the luminosity distance D_L (z) as D(z) = D_L (z)/(1 + z) (Hogg 1999), and ${D}^{{\prime} }(z)={dD}(z)/{dz}$ denotes the derivative of D(z) with respect to redshift z. More specifically, we use the Gaussian processes (GP) method (Seikel et al. 2012) to reconstruct the first derivative of luminosity distance D_L (z). This method, which assumes that the distribution of data is Gaussian, can effectively reconstruct a function and its derivatives from a given data set without parameterization (Shafieloo et al. 2012; Cao et al. 2019a; Qi et al. 2019a; Liu et al. 2019; Wu et al. 2020; Zheng et al. 2020). The reconstruction of ${D}_{L}^{\prime} (z)$, together with the observations of D_L (z) and X(z) at individual redshifts, will provide different measurements of Ω_k through Equation (16).

It is necessary to give a brief introduction to the Gaussian process, which is executed by the Python package Gapp ¹⁰ in this work. Given a data set ${ \mathcal D }$ of observations: ${ \mathcal D }=\{({x}_{i},{y}_{i})\}$ (in the cases investigated in this study, x and y are the redshift z and the luminosity distance D_L from the simulation, respectively). The covariance function cov(f(x), f($\tilde{x}$)) = k(x, $\tilde{x}$) is used to describe the connection between the function value at x and the function value at another point $\tilde{x}$. Here, we consider the squared exponential covariance function

$$\begin{eqnarray}&&k(x,\tilde{x})={\sigma }_{f}^{2}\exp \left(-\displaystyle \frac{{\left(x-\tilde{x}\right)}^{2}}{2{{\ell }}^{2}}\right),\end{eqnarray} \tag{ 17 }$$

where the characteristic lengths ℓ and σ_f represent the typical changes in f(x) in the x-direction and y-direction, respectively. These two hyperparameters can be trained by the observational data.

First, we reconstruct a function f(x) to describe the data set, which can be expressed by the mean value μ(x) and the covariance function cov(f(x), f($\tilde{x}$)):

$$\begin{eqnarray}&&f(x)\sim { \mathcal G }{ \mathcal P }\left(\mu (x),k(x,\tilde{x})\right).\end{eqnarray} \tag{ 18 }$$

However, there is a difference between the real data y and the reconstructed Gaussian function f(x): y_i = f(x_i ) + _i , where Gaussian noise _i with variance ${\sigma }_{i}^{2}$ is assumed. Therefore, for a set of observational points $X=\left\{{x}_{i}\right\}$, the observational data can be written as

$$\begin{eqnarray}&&y\sim { \mathcal G }{ \mathcal P }(\mu ,K(X,X)+C),\end{eqnarray} \tag{ 19 }$$

by adding the variance to the covariance matrix C of the data.

For observations with a limited sample size, one may need to reconstruct another function f* to describe the observational data but at some other points X*, which are typically an extended point set. Simply, we obtain

$$\begin{eqnarray}&&f* \sim { \mathcal G }{ \mathcal P }(\mu * ,K(X* ,X* )),\end{eqnarray} \tag{ 20 }$$

where μ* is a prior assumed mean of f*. Combining these two Gaussian processes for y and f* (Equations (19) and (20)) and calculating the conditional distribution, we can reconstruct the mean and covariance of f* by

$$\begin{eqnarray}&&\overline{f* }=\mu * +K(X* ,X){\left[K(X,X)+C\right]}^{-1}(y-\mu )\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{cov}(f* )=K(X* ,X* )-K(X* ,X)\\ \,\,\times \,{\left[K(X,X)+C\right]}^{-1}K(X,X* ).\end{array}\end{eqnarray} \tag{ 22 }$$

And the variance of f* can be simply obtained by diagonalizing cov(f*). Thus, we can expand the observational data set.

The derivative of the function f* can also be calculated by giving the Gaussian process for y and $f* ^{\prime}$. Similarly, we have

$$\begin{eqnarray}&&\overline{f* ^{\prime} }=\mu * ^{\prime} +K^{\prime} (X* ,X){\left[K(X,X)+C\right]}^{-1}(y-\mu )\end{eqnarray} \tag{ 23 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{cov}(f* ^{\prime} ) & = & -K^{\prime} (X* ,X){\left[K(X,X)+C\right]}^{-1}K^{\prime} (X,X* )\\ & & +K^{\prime\prime} (X* ,X* ),\end{array}\end{eqnarray} \tag{ 24 }$$

where

$$\begin{eqnarray}&&{\left[K^{\prime} (X,X* )\right]}_{{ij}}=\mathrm{cov}\left({f}_{i},\displaystyle \frac{\partial f{* }_{j}}{\partial x{* }_{j}}\right)=\displaystyle \frac{\partial k({x}_{i},x{* }_{j})}{\partial x{* }_{j}}\end{eqnarray} \tag{ 25 }$$

is the covariance between the function f* and its derivative, and

$$\begin{eqnarray}&&{\left[K^{\prime\prime} (X* ,X* )\right]}_{{ij}}=\mathrm{cov}\left(\displaystyle \frac{\partial f{* }_{i}}{\partial x{* }_{i}},\displaystyle \frac{\partial f{* }_{j}}{\partial x{* }_{j}}\right)=\displaystyle \frac{{\partial }^{2}k(x{* }_{i},x{* }_{j})}{\partial x{* }_{i}\,\partial x{* }_{j}}.\end{eqnarray} \tag{ 26 }$$

is the covariance between its derivatives. Then both function f* and its derivative $f* ^{\prime}$ can be reconstructed from the observational data ${ \mathcal D }$.

We first combine D_L (z) and X(z) to provide a numerical constraint on Ω_k by calculating the χ² statistics in Equation (15). The result from 10,000 simulated GW events detected by DECIGO is

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{k}=-0.05\pm 0.12.\end{eqnarray} \tag{ 27 }$$

This result is consistent with the fiducial value of Ω_k = 0 within a 1σ confidence level. Compared to the other model-independent tests involving H(z) and distances from popular cosmological probes, our result is more precise than that of Ω_k = 0.63 ± 0.34 from the Pantheon SNe Ia plus H(z) (CC) (Wang et al. 2020), and the result Ω_k = −0.92 ± 0.43 from UV + X-ray quasars plus H(z) (CC) (Wei & Melia 2020).

In order to demonstrate the precision of the assessment of curvature parameter with a certain number of GW events, we show the best-fit Ω_k and 1σ confidence level as a function of the number of GW events in Figure 3 and Table 1. The model-independent test of Ω_k from Pantheon SNe Ia (blue diamond), which will be introduced in the later analysis, is also plotted for comparison. As one may see, the precision of the determined Ω_k from 6000 GW events (ΔΩ_k = 0.165) is more competitive than that from the Pantheon SNe Ia sample.

Zoom In Zoom Out Reset image size

Table 1. Summary of Model-independent Curvature Determinations from Different Numbers of GW Events

N	Ωk	N	Ωk
1000	${0.03}_{-0.37}^{+0.52}$	6000	−0.03 ± 0.17
2000	0.01 ± 0.29	7000	−0.07 ± 0.15
3000	−0.05 ± 0.22	8000	−0.05 ± 0.14
4000	−0.03 ± 0.20	9000	−0.01 ± 0.13
5000	−0.05 ± 0.18	10,000	−0.05 ± 0.12

Download table as:  ASCIITypeset image

The individual measurements of 10,000 Ω_k (z) obtained from Equation (16) in Section 2.4 are shown in Figure 4. To estimate the improvement of the combined measurement, in Figure 4 we also summarize the multiple Ω_k within the redshift bin Δz = 0.1 through inverse variance weighting, which allows a direct check of its predicted constancy with the given redshift. We find that with increasing redshift the derived Ω_k remains within the error bar (68.3% confidence level) of the flat case, which underlies the assumption of our GW data simulations. It is worth noticing that the uncertainty of Ω_k , which is much higher in the low redshift range, also fluctuates at high redshifts (3 < z < 5). Such a tendency could be explained by the term ${\left({D}^{{\rm{{\prime} }}}(z)/D(z)\right)}^{2}E(z){\sigma }_{E(z)}$ in the Ω_k (z) error equation derived from Equation (16), which dominates the uncertainty of Ω_k measurements. In this term, the function ${\left(D^{\prime} (z)/D(z)\right)}^{2}$ with large values at z ∼ 0 tends to decrease with increasing redshift, which generates a relatively small Ω_k uncertainty at higher redshifts. However, the function E(z)σ_E(z) exhibits the opposite tendency, generating a fluctuation of Ω_k uncertainty at higher redshifts when combined with the function ${\left(D^{\prime} (z)/D(z)\right)}^{2}$.

Zoom In Zoom Out Reset image size

For comparison, we also provide model-independent measurements of Ω_k from electromagnetic observations. The 1048 luminosity distance measurements from Pantheon SNe Ia, as well as 41 OHD from CC and BAO measurements (Gaztañaga et al. 2009; Blake et al. 2012; Busca et al. 2013; Samushia et al. 2013; Xu et al. 2013; Font-Ribera et al. 2014; Delubac et al. 2015), are used to test the cosmic curvature. Let us note that in order to achieve a better redshift match with the supernova data, we use the GP to obtain the reconstruction of a smooth H(z) function, with the generation of 200 reconstructed H(z) data well matched to the redshifts of the Pantheon SNe Ia sample.

According to the BEAMS with Bias Corrections (BBC) method (Kessler & Scolnic 2017), the observed distance modulus of SNe can be simply given by the apparent magnitude (m_B ) and absolute B-band magnitude (M_B ) as μ = m_B − M_B , with the nuisance parameters in the Tripp formula (Tripp 1998) retrieved. The theoretical distance modulus m_th can then be obtained from

$$\begin{eqnarray}&&{m}_{\mathrm{th}}=5\mathrm{log}\displaystyle \frac{{D}_{L}^{\mathrm{th}}}{\mathrm{Mpc}}+25+{M}_{B},\end{eqnarray} \tag{ 28 }$$

where ${D}_{L}^{\mathrm{th}}$ is given by Equation (2) involving the reconstructed H(z). Then, the cosmic curvature Ω_k and the parameter M_B can be constrained by minimizing the χ² statistic

$$\begin{eqnarray}&&{\chi }_{\mathrm{SNe}}^{2}=\sum _{i=1}\displaystyle \frac{{\left[{m}_{i}^{\mathrm{obs}}-{m}_{i}^{\mathrm{th}}\right]}^{2}}{{\sigma }_{\mathrm{SNe}}^{2}+{\sigma }_{{m}_{\mathrm{th}}}^{2}},\end{eqnarray} \tag{ 29 }$$

where σ_SNe accounts for the error in SNe Ia observations, propagated from the covariance matrix in Scolnic et al. (2018). The marginalized probability distribution of each parameter and the marginalized 2D confidence contours are shown in Figure 5. The best-fit cosmic curvature and the absolute B-band magnitude with 1σ are Ω_k = −0.16 ± 0.17 and M_B = −19.32 ± 0.01. This result favors a zero cosmic curvature at 68.3% confidence level, which shows no evidence for a deviation from a flat universe at the current observational data level. Meanwhile, the uncertainty of Ω_k from the GW method (ΔΩ_k ∼ 0.12) is ∼30% smaller than that from the SNe Ia method. Moreover, the SNe Ia method is strongly dependent on the choice of the reconstruction methods of the H(z) function, as was pointed out in the recent work of Wang et al. (2020). The determined cosmic curvature Ω_k is strongly degenerate with the absolute magnitude M_B of SNe Ia, similar to the results obtained by examining the cosmic opacity with gravitational waves and SNe Ia (Qi et al. 2019c).

Zoom In Zoom Out Reset image size

Then, we transform the distance modulus μ to D_L through

$$\begin{eqnarray}&&{D}_{L}(z)={10}^{\mu (z)/5-5}\,({\rm{M}}{\rm{p}}{\rm{c}}),\end{eqnarray} \tag{ 30 }$$

where the absolute magnitude of SNe Ia is set at the best-fit value M_B = −19.32 from the results shown in Figure 5. Following the same procedure as the GW observations, the cosmic curvature is measured at 41 different redshifts, as shown in Figure 6. Compared with the GW method, the redshift of Ω_k is only up to z ∼ 2.3, limited by the redshift coverage of SNe Ia (0.01 < z < 2.26) and OHD (0.07 < z < 2.36). Meanwhile, the Ω_k measurement based on our GW method is more precise than those from the current EM observations, which indicates another advantage of our methodology in testing the spatial curvature in the GW domain.

Zoom In Zoom Out Reset image size

The performance of the current traditional EM method is restricted by the big gap between the sample sizes, as well as the redshift coverage of the observational H(z) data and D_L (z) data. On the one hand, the Hubble diagram of SNe Ia contains only of the order of ∼10³ SNe Ia. Such a situation will be greatly improved in the era of the LSST, which could discover an unprecedented number of SNe Ia (∼10⁶), with a large fraction (∼10%) expected to be turned into useful distance indicators (Lochner 2022). On the other hand, future observations of redshift drift (Sandage 1962), which is also known as the Sandage–Loeb test, provide an important method to derive precise measurements of H(z) at different redshifts. Specifically, by observing the redshift drift in the optical and radio bands, the E-ELT and the Square Kilometre Array will offer H(z) measurements in the redshift ranges 2 < z < 5 and 0 < z < 0.3, respectively (Liske et al. 2008; Quercellini et al. 2012; Martins et al. 2016). In addition, strongly lensed SNe Ia, which will also be discovered in larger numbers by LSST, enable a more precise model-independent probe of cosmological parameters based on the distance sum rule (Cao et al. 2018; Ma et al. 2019b). In particular, based on the simulated sample of 200 lensed SNe Ia with time-delay measurements (Qi et al. 2022), model-independent constraints on the Hubble constant H₀ and cosmic curvature parameter Ω_k would be achieved with high precision (ΔH₀ = 0.33 km s⁻¹ Mpc⁻¹ and ΔΩ_k = 0.05). Therefore, one might be optimistic about achieving much higher precision of improved EM observations in the next decades. In that case, the precision of Ω_k measurements in the EM domain would be much higher. On the other hand, considering the difficulty of deriving multiple measurements (D_L (z) and H(z)) at exactly the same redshift, the prospects for constraining the cosmic curvature in the GW domain could be promising, based on the combination of two different observables for the same objects at high redshifts.

As a standard siren, GWs from a binary neutron star merger provide a direct way to measure the luminosity distance (D_L ) without the need for a cosmological distance ladder. In addition, the accelerating expansion of the universe may cause an additional phase shift in the gravitational waveform, which would allow us to measure the acceleration parameter. Thus, GW measurement provides an important opportunity to determine the curvature parameter Ω_k in the GW domain based on the combination of two different observables for the same object at high redshifts.

In this paper, we investigate how such an idea could be implemented with the future generation of the space-based DECIGO in the framework of two model-independent methods. Our results show that DECIGO could provide a reliable and stringent constraint on the cosmic curvature at a precision of ΔΩ_k = 0.12, which is comparable to the latest model-independent estimations using different EM probes. Furthermore, we use the GP method to reconstruct the first derivative of D_L and obtain 10,000 individual measurements of Ω_k at different redshifts (z ∼ 5). Compared to the traditional model-independent estimations of the spatial curvature using other EM observations, GW sirens have several benefits as follows.

• 1.  
  Cosmological-model-independent: measurements of GW standard sirens could provide independent measurements of luminosity distance and acceleration parameter without the cosmological distance ladder. And the matched filtering analysis in the GW method does not require the assumption of any fiducial cosmological model. In addition, the Ω_k from EM method is strongly degenerate with other parameters such as the absolute magnitude, while that from GW standard sirens has the advantage of no nuisance parameters.
• 2.  
  High redshift: GW detectors can observe a large number of events at high redshifts, which allows us to probe the cosmic curvature at high redshifts. However, the current OHD obtained by radial BAOs and CCs still provide no high-redshift data.
• 3.  
  Well redshift-matched: both the luminosity distance and the acceleration parameter can be determined from each GW event, which means that they are already well matched in redshift and can be used directly for the statistical constraint on Ω_k (as in Equation (15)). In the traditional EM method, due to the large difference in the amount of observational H(z) data and observational D_L (z) data, a parametric or nonparametric method must be used to reconstruct H(z) first, restricting the constraint ability of SNe Ia. In addition, the redshift match between two sets of data also leads to errors.

Summarizing, the GW observations provide a powerful and novel method to estimate the spatial curvature in different cosmological-model-independent ways (Zheng et al. 2021; Cao et al. 2022b). This strengthens the probative power of our method, especially in the framework of DECIGO, to inspire other new observing programs and theoretical works in the near future. However, there are still some issues that should be emphasized here.

• 1.  
  The redshifts of the GW sources are necessary ingredients to derive model-independent constraints on the cosmic curvature with DECIGO. Considering the angular resolution of DECIGO (∼1 arcsec²), which can uniquely identify the host galaxy of the binary, we could adopt a widely used method such as the optical (or infrared) identification of the host galaxy of the GW event. For high-redshift GWs, one could turn to high-redshift tracers such as quasars or gamma-ray bursts (GRBs), along with follow-up observations to determine their redshifts, considering the significant development of multimessenger astronomy in the framework of DECIGO. In addition, several methods have been proposed in the literatures to address this issue, such as the "galaxy voting" method (redshift distribution for host galaxies) (MacLeod & Hogan 2008; Trott & Huterer 2021), the redshift distribution for coalescing sources (Ding et al. 2019), neutron star mass distribution (Taylor et al. 2012; Taylor & Gair 2012), cross-correlation of gravitational-wave standard sirens and galaxies (Oguri 2016; Mukherjee et al. 2021, 2022) and the tidal deformation of neutron stars (Messenger & Read 2012; Messenger et al. 2014; Wang et al. 2020).
• 2.  
  Different from current observational H(z) data, the measurement of acceleration parameter X is based on time measurement in the observer coordinate (which is similar to the measurement of redshift drift), while the OHD rely on time measurement in the universe coordinate. And X appears in the 4PN order correction in the GW waveform, which requires high-precision detection of a GW signal especially at lower frequencies (Seto et al. 2001; Nishizawa et al. 2012). Fortunately, as one of its major objectives, DECIGO is capable of direct measurement of the acceleration of the universe. In the future, DECIGO will detect a large number of NS binaries in inspiraling phases, which will provide an unprecedented opportunity for high-precision detections of cosmic acceleration and will open up a window for gravitational-wave cosmology.

This work was supported by the National Natural Science Foundation of China under grant Nos. 12021003, 11690023, and 11920101003; the National Key R&D Program of China No. 2017YFA0402600; the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23000000; and the Interdiscipline Research Funds of Beijing Normal University.