Multiple measurements of gravitational waves acting as standard probes: model-independent constraints on the cosmic curvature with DECIGO

Although the spatial curvature has been precisely determined via the cosmic microwave background (CMB) observation by Planck satellite, it still suffers from the well-known cosmic curvature tension. As a standard siren, gravitational waves (GWs) from binary neutron star mergers provide a direct way to measure the luminosity distance. In addition, the accelerating expansion of the universe may cause an additional phase shift in the gravitational waveform, which allows us to measure the acceleration parameter. This measurement provides an important opportunity to determine the curvature parameter $\Omega_k$ in the GW domain based on the combination of two different observables for the same objects at high redshifts. In this study, we investigate how such an idea could be implemented with future generation of space-based DECi-hertz Interferometer Gravitational-wave Observatory (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 $\Delta\Omega_k$=0.12, which is comparable to existing results based on different electromagnetic data. Our constraints are more stringent than the traditional electromagnetic method from the Pantheon SNe Ia sample, which shows no evidence for the deviation from the flat universe at $z\sim 2.3$. More importantly, with our model-independent method, such a second-generation space-based GW detector would also be able to explore the possible evolution $\Omega_k$ with redshifts, through direct measurements of cosmic curvature at different redshifts ($z\sim 5$). Such a model-independent $\Omega_k$ reconstruction to the distance past can become a milestone in gravitational-wave cosmology.


INTRODUCTION
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 Friedmann-Lemaître-Robertson-Walker (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 closed to zero (Weinberg et al. 2013). Current cosmological observations, particularly the latest Planck2018 results, which combined the cosmic microwave background (CMB) and baryon acoustic oscillation (BAO) measurements, strongly favor a flat universe: Ω k = 0.0007±0.0019 (Planck Collaboration, et al. 2018) (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: Ω k = −0.044 +0.018 −0.015 (Planck Collaboration, et al. 2018;Di Valentino, Melchiorri & Silk 2020). In addition, such a stringent constraint on the cosmic curvature strongly relies 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 model Clarkson, Cortes & Bassett 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 modelindependent ways.
The distance sum rule has been proposed as a model-independent method to constrain the curvature of the universe, which is generally implemented using strong lensing observations (Cao & Zhu 2012a;Cao et al. 2012bCao et al. , 2015Ma et al. 2019a) with other distance measurements, such as type Ia supernovae (SNe Ia) (Räsänen, Bolejko & Inogueno 2015;Denissenya, Linder & Shafieloo 2018) or quasar (QSO) (Qi et al. 2019b;Zhou et al. 2020;Lian et al. 2021). Another model-independent way to constrain the cosmic curvature is 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, Bassett & Lu 2008;Shafieloo & Clarkson 2010;Mortsell & Jonsson 2011;Sapone, Majerotto & Nesseris 2014;Cai, Guo & Yanget 2016). This test has been fully implemented with various observational data. For example,  recently used the latest Pantheon SNe Ia data combined with cosmic chronometers (CC) to constrain the cosmic curvature; a well-measured quasar sample (Risaliti & Lusso 2015;Lusso & Risaliti 2016;) also shows a great potential for probing the cosmic curvature (Wei et al. 2020);Cao et al. (2019b) proposed an improved model-independent test of cosmic curvature with ultra-compact structures in radio quasars as standard rulers, which shows which shows no evidence for the deviation from the flat universe; and Takada & DoréO. (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, Cortes & Bassett 2007), which also requires the observational data of 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 spacial 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 catalogue 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), provide an absolute measurement of the luminosity distance. Only if these binary mergers are accompanied by short-hard γ-ray bursts (shGRB), can they be observed through both electromagnetic (EM) and GW. 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 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 gravitational-wave signals from a neutron star binary at a redshift of 1, the universe expansion acceleration would cause a phase delay of only 1 sec during the 10-year probe (Seto et al. 2001;Nishizawa 2012). As a future space-borne GW detector, DECIGO (Kawamura et al. 2011;Seto et al. 2001) is designed to improve the detection sensitivity of GW at lower frequencies, with its most sensitive frequencies between 0.1 and 10 Hz. DECIGO will have four clusters of spacecraft, and each cluster consists 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 −25 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 GW 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 ranges of 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 observational Hubble parameter data (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.

METHODOLOGY
Assuming that the universe is homogeneous and isotropic on large scales, it can be described by the FLRW metric: 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 matter and energy of the universe. The dimensionless curvature K = −1, 0, +1 corresponds to open, flat and closed universes respectively. With such a metric, the luminosity distance D L (z) can be expressed as for Ω k > 0, for Ω K = 0, for Ω k < 0.
(2) The dimensionless Hubble parameter E(z) is defined as H(z)/H 0 , where H(z) is the universe expansion rate and H 0 is the Hubble constant. The curvature parameter Ω k is related to K as Ω k = −c 2 K/(a 0 H 0 ) 2 , where c is the speed of light.
2.1. Measuring the distance and the acceleration parameter with GWs As a standard siren, 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 neutronstar 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 1 and m 2 . One can define the chrip mass M c = M η 3/5 , and the redshifted chirp mass M z ≡ M (1 + z c )η 3/5 , where M = m 1 + m 2 is the total mass of the binary system, η = m 1 m 2 /M 2 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 ( 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 1 : 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 2 By substituting Eq. (4) into Eq.
(3) and applying the stationary phase approximation (Cutler & Flanagan 1994), we can transform the Fourier component of the waveform into the frequency domaiñ where corresponds to the gravitational waveform with cosmic acceleration, with x ≡ (πM z f ) 2/3 and Ψ N (f ) ≡ 3 128 (πM z f ) −5/3 . The gravitational waveform without cosmic acceleration is given bỹ with the amplitude written as 1 Such a relation is only an approximation, without considering the contribution from high-order terms.
2 Note that X(z) is related to the redshift drift ∆tz as ∆tz = The polarisation amplitude A pol,α (t), the polarisation phases ϕ pol,α (t) (α = I, II represents the number of individual detectors 3 ) and the Doppler phase ϕ D (t) are given in Yagi & Tanaka (2010). For the phase of Ψ(f ), we use the restricted-2PN waveform including spin-orbit coupling at the 1.5PN order 4 (Maggiore 2008;Kidder et al. 1993), which is also reported by in 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 Eq. (6), we take the binary parameters as where β is the spin-orbit coupling parameter, and φ c represents the coalescence phase. (θ S ,φ S ) indicates the direction of the source in the solar barycenter frame, and (θ L ,φ L ) specifies the direction of the orbital angular momentum.
One can use Fisher analysis to estimate the measurement accuracies of the binary parameters θ i . The measurement accuracy is given by ∆θ i ≡ (Γ −1 ) ii (Cutler & Flanagan 1994), with the Fisher matrix: where S h (f ) is the analytical expression of the DE-CIGO noise power spectrum (Kawamura et al. 2006(Kawamura et al. , 2019Yagi & Seto 2011). The lower cut-off of frequency corresponds to the frequency at which coalescence begins to be observed, with ∆T obs representing the observation time. And f max is the higher cut-off frequency of the detector and is set equal to 100 Hz. In addition, for DE-CIGO, there will be eight uncorrelated interferometric signals (Kawamura et al. 2011;Yagi & Seto 2011); thus, the Fisher matrix above should be multiplied by 8.
For the convenience of these calculations, we set m 1 = m 2 = 1.4M ⊙ and take t c = φ c = β = 0. For each fiducial redshift z c , we randomly generate 10 4 sets of (θ S ,φ S ,θ L ,φ L ), and for each set, we calculate the uncertainty (Γ −1 ) 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 2-Dimension submatrix of D L and X, with the instrumental uncertainty of X being: Similarly, the measurement error of luminosity distance is estimated as Considering the lensing uncertainty caused by the weak lensing effect σ lens DL = 0.05zD L (Sathyaprakash et al.  2010), the luminosity distance error is considered to be: For the redshifts of GW sources, we sample them from the merger rate of double compact objects that reflects the star formation history (Dominik et al. 2013), taking the data from the so-called "rest frame rates" in the cosmological scenario 5 . In our simulation, we assume a flat ΛCDM as the proposed fiducial model with the cosmological parameters derived by Planck2018 measurements (Planck Collaboration, et al. 2018). The central values of the simulated 10,000 X(z) and D L (z) measurements are shown in Fig. 1, and the acceleration parameters and their corresponding errors are shown in Fig. 2 Our main concern here is what fraction of binary sources can we determine the redshifts in the era of DECIGO. For neutron-star (NS) binaries, the most reliable method is determining the electromagnetic counterpart of GW events directly. Meanwhile, benefit from the angular resolution of ∼ 1 arcsec 2 , DECIGO is expected to uniquely identify the host galaxy of the binary (Cutler & Holz 2009), the redshift of which could be determined from multi-messenger EM observations.
We now estimate the number of potential host galaxies with redshift determination for the binaries. Following the methodology proposed in Holz & Hughes (2005), there would be more than ∼ 10 11 galaxies potentially observable over the entire celestial sphere, with the number density of galaxies ∼ 10 3 /arcmin 2 inferred from the observation of the Hubble Ultra Deep Field (HUDF) (Beckwith et al. 2006;Cutler & Holz 2009). While the future galaxy spectroscopic surveys, JDEM/WFIRST (Gehrels 2010;Michael et al. 2011), plan to obtain spectroscopic redshifts for 10 8 galaxies in the redshift range 0.5< z <2 with a precision better than 0.1% (the number of galaxies with photometric redshift measurements is expected to be larger, ∼ 10 9 ). Thus, the fraction of galaxies whose redshifts are listed in the galaxy catalogue is ∼ 10 −3 . DECIGO is expected to observe 10 6 GW events coming from the neutron star binaries within redshift z ∼ 5 (Kawamura et al. 2019;Nishizawa 2012). If the GW events randomly occur in any one of the galaxies, there would be ∼ 10 3 events with redshift determination (∼ 10 4 if photometric galaxies are considered). In addition, the proposed wide-field survey BigBOSS (Schlegel et al. 2010) would measure ∼ 5 million spectroscopic redshifts per year for galaxies in the range 0.2< z <3.5. In the era of DECIGO, LSST may have already determined the photometric redshifts for a large fraction of host galaxies in ∼ 1/3 of the sky (Cutler & Holz 2009). Such methodology could be easily extended to z ∼ 3, focusing on Euclid's near-infrared photometry combined with ground-based optical photometry.
For high-redshift GWs, one could turn to high-redshift tracers such as quasars or gamma-ray bursts (GRBs). The Gamow Explorer program, which has been proposed for searching for X-ray and optical-IR counterparts to high-redshift GW events (White et al. 2021), could rapidly detect the GRB (by Lobster Eye X-ray Telescope) and provide GRB redshift estimates (by Photo-z Infra Red Telescope). It also enables space and groundbased observatories to follow up this GRB to determine the redshift of its host galaxy and study the afterglow in detail 6 . For the follow-up observations, the deep field of HST is useful to observe galaxies that are fainter than the characteristic brightness, which not only contributes to the observations of galaxies at a very high redshift (z ∼ 11), but also provides important observations for galaxies in the redshift range of 2< z <5 (Beckwith et al. 2006). As a scientific successor to HST, James Webb Space Telescope (JWST) (Gardner et al. 2006) can detect galaxies at much higher redshifts, and also observe faint infrared afterglows of short GRBs from BNS merg-ers at a distance of 150 Mpc (Lu et al. 2021), and kilonovae within a distance of ∼200 Mpc (Bartos et al. 2016). Other ground-based facilities would also be alerted a few years before the mergers, including Keck Observatory, Gran Telescopio Canarias (GTC), Gemini, Very Large Telescope array (VLT) (Hartoog et al. 2015) and future planned 40 m facilities such as E-ELT. We remark here that DECIGO can observe GW signals several years before the coalescence, and predict coalescence time with an accuracy of ∼ 0.1 sec (Kawamura et al. 2021). Based on the precise time and sky location of GW events months in advance, simultaneous gamma-ray observations and electromagnetic follow-up observations would be more reliable and frequent. Therefore, multi-messenger astronomy will develop significantly under the DECIGO framework (Chen et al. 2021;Cao et al. 2022).
In the above estimation, we only provide the fraction of galaxies whose redshifts are listed in the galaxy catalogue. Such a worst case with large uncertainty could be markedly improved, with dedicated follow-up observations targeted at the GW events or the host galaxies. Actually, the fraction of neutron-star binaries with redshift determination would be much larger, considering the fact that GW events are likely to occur in more massive luminous galaxies that are easier to observe. In this sense, it is reasonable to consider a conservative case during analysis, with the simulation of 10,000 GW events used to derive numerical constraints and redshift reconstruction of cosmic curvature.
2.3. Numerical constraints on cosmic curvature Ω k (z) GW 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 X L (z) inferred from Eq. (2), where curvature parameter Ω k is considered and E(z) can be derived from the definition equation of X(z) in Eq. (5), to provide a numerical constraint on the cosmic curvature Ω k (z). This process could be achieved by minimizing the χ 2 statistic where σ D X L can be derived from Eq. (2) and σ DL is given in Eq. (14). We use the Markov Chain Monte Carlo (MCMC) method (Foreman-Mackey & Hogg 2013) to obtain the best-fit value of Ω k and its uncertainty, which is a model-independent estimation of the cosmic curvature.

Individual measurements of cosmic curvature Ω k (z)
Future space-based GW detectors, such as DECIGO, will aim to detect a large number of neutron-star binaries at much higher redshifts, which enables the reconstruction of cosmic curvature in the early universe. The cosmic curvature Ω k can be expressed as (Clarkson, Cortes & Bassett 2007) where the expansion rate of universe H(z) can be obtained from the GW measurement of X(z), the trans-verse comoving distance D(z) is simply related to the luminosity distance D L (z) as D(z) = D L (z)/(1 + z) (Hogg et al. 1999), and D ′ (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 derivation 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 parametrization (Shafieloo, Kim & Linder 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 (z), together with the observations of D L (z) and X(z) at individual redshifts, will provide different measurements of Ω k through Eq. (16).
It is necessary to give a brief introduce to the Gaussian process, which is executed by the Python package Gapp 7 in this work. Given a data set D of observations: 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 (x)) = k(x,x) is used to describe the connection between the function value at x and the function value at the other pointx. Here, we consider the squared exponential covariance function where the characteristic lengths ℓ and σ f represent the typical changes in f (x) in 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 (x)): 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 σ 2 i is assumed. Therefore, for a set of observational points X = {x i }, the observational data can be written as 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 where µ * is a prior assumed mean of f * . Combining these two Gaussian processes for y and f * (Eq. (19) and Eq. (20)) and calculating the conditional distribution, we can reconstruct the mean and covariance of f * by and cov(f * ) = K(X * , X * )−K(X * , X)[K(X, X)+C] −1 K(X, X * ).
(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 * ′ . Similarly, we have and where is the covariance between the function f * and its derivative, and (26) is the covariance between its derivatives. Then both function f * and its derivative f * ′ can be reconstructed from the observational data D.

RESULTS AND DISCUSSION
We first combine the D L (z) and X(z) to provide a numerical constraint on Ω k by calculating the χ 2 statistics in Eq. (15). The result from 10,000 simulated GW events detected by DECIGO is This result is consistent with the fiducial value of Ω k = 0 within a 1σ confidence level. Compared to the other model-independent tests involving the H(z) and distances from popular cosmological probes, our result is more precise than that of the Pantheon SNe Ia plus H(z) (CC) as Ω k = 0.63 ± 0.34 , and the result of UV+Xray quasars plus H(z) (CC) as Ω k = −0.92 ± 0.43 (Wei et al. 2020). In order to demonstrate the precision of the curvature parameter assessment 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 Fig. 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 of Pantheon SNe Ia sample.
The individual measurements of 10,000 Ω k (z) obtained from Eq. (16) in Section 2.4 are shown in Fig. 4. To estimate how much improvement of the combing measurement, in Fig. 4 we also summarize the multiple Ω k within the redshift bin of ∆z = 0.1 through the inverse variance weighting, which allows a direct check  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 10000 −0.05 ± 0.12 of its predicted constancy with the given redshift. We find that with the increase of redshift, the derived Ω k remains within the error bar (68.3% confidence level [C. L.]) 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 [(D ′ (z)/D(z)) 2 E(z)σ E(z) ] in the Ω k (z) error equation derived from Eq. (16), which dominates the uncertainty of Ω k measurements. In this term, the function of (D ′ (z)/D(z)) 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 of [E(z)σ E(z) ] exhibits an opposite tendency, generating a fluctuation of Ω k uncertainty at higher redshifts when combined with the function of (D ′ (z)/D(z)) 2 .
For comparison, we also provide model-independent measurements of Ω k from electromagnetic observations. The 1048 luminosity distance measurements from SNe Ia Pantheon, as well as 41 observational Hubble parameter data (OHD) from cosmic chronometer (CC) and baryon acoustic oscillation (BAO) measurements (Gaztanaga, Cabré & Hui 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 GP to obtain the reconstruction of a smooth H(z) function, with the generation of 200 reconstructed H(z) data well matched the redshifts of Pantheon SNe Ia sam-

ple.
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 (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 by where D th L is given by Eq.
(2) involving the reconstructed H(z). Then, the cosmic curvature Ω k and the parameter M B can be constrained by minimizing the χ 2 statistic where σ SN e 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 2-D confidence contours are shown in Fig. 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 the deviation from the 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 de- pendent on the choice of the reconstruction methods of H(z) function, as was pointed out in the recent works of . The determined cosmic curvature Ω k is strongly degenerated with the absolute magnitude M B of SNe Ia, similar to the results obtained by examining the cosmic opacity with gravitational waves and Type Ia Supernova (Qi et al. 2019c). Then, we transform the distance modulus µ to D L through D L (z) = 10 µ(z)/5−5 (Mpc), where the absolute magnitude of SNe Ia is set at the bestfit value M B = −19.32 from the results shown in Fig. 5. Following the same procedure as the GW observations, the cosmic curvature is measured at 41 different redshifts, as shown in Fig. 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.
The performance of the current traditional EM method is restricted by the big gap between the amount, 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 type Ia supernovae (SNe Ia) contains only the order of ∼ 10 3 SNe Ia. Such situation will be greatly improved in the era of the Large Synoptic Survey Telescope (LSST), which could discover an unprecedented number of SNe Ia ∼ 10 6 , with a large fraction (∼ 10%) expected to be turned into useful distance indicators (Lochner et al. 2021). 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 European Extremely Large Telescope (E-ELT) and the Square Kilometre Array (SKA) will offer the H(z) measurements in the redshift ranges of 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, enables a more precise model-independent probe of cosmological parameters based on the distance sum rule (Cao et al. 2018;Ma et al. 2019b). Specially, based on the simulated sample of 200 lensed SNe Ia with time-delay measurements (Qi et al. 2022), model-independent constraints of the Hubble constant H 0 and cosmic curvature parameter Ω k would be achieved with high precision (∆H 0 = 0.33 kms −1 M pc −1 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. Whereas, considering the the difficulty of deriving multiple measurements (D L (z) and H(z)) exactly at 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.

CONCLUSIONS
As a standard siren, gravitational wave (GW) from binary neutron star merger provides a direct way to measure the luminosity distance (D L ) without the need of cosmological distance ladder. In addition, the accelerating expansion of the universe may cause an additional phase shift in the gravitational waveform, which allows 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 objects at high redshifts.
In this paper, we investigate how such an idea could be implemented with future generation of space-based DECi-hertz Interferometer Gravitational-wave Observatory (DECIGO) in the framework of two modelindependent 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 ob-tain individual 10,000 measurements of Ω k at different redshifts (z ∼ 5). Compared to the traditional modelindependent estimations of the spatial curvature using other EM observations, GW sirens have several benefits as follows: • Cosmological-model independent: GW standard siren measurements could provide independent measurements of luminosity distance and acceleration parameter without the cosmological distance ladder. And the matched filtering analysis in GW method does not require an 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.
• 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 baryon acoustic oscillations (BAO) and cosmic chronometer (CC) still provide no high redshift data.
• 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 Eq. (15)). While the traditional EM method, due to the large gap between the amount of observational H(z) data and observational D L (z) data, a parametric or non-parametric method must be used to reconstruct H(z) before, restricting the constraint ability of SNe Ia. In addition, the redshift match between two sets of data also leads to errors.
However, there are still some issues that should be emphasized here: • The redshifts of the GW sources are necessary ingredients to perform model-independent constraints on the cosmic curvature with DECIGO.
Considering the angular resolution of DECIGO (∼ 1 arcsec 2 ), which can uniquely identify the host galaxy of the binary, we could take a widely used method that through the optical (or infrared) identification of the host galaxy of the GW event.
For high-redshift GWs, one could turn to highredshift tracers such as quasars or gamma-ray bursts (GRBs), along with follow-up observations to determine their redshifts, considering the significant development of multi-messenger 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 et al. 2008;Trott & Huterer 2021), the redshift distribution for coalescing sources (Ding et al. 2021), neutron star mass distribution (Taylor et al. 2012a,b), cross-correlation of gravitational wave standard sirens and galaxies (Oguri et al. 2016), and the tidal deforma-tion of neutron stars (Messenger & Read 2012;Messenger & Takami 2014;Wang et la. 2020).
• Different from current observational H(z) data, the measurement of acceleration parameter X is based on time measurement in observer coordinate (which is similar to the measurement of redshift drift), while the OHD relies on the time measurement in the university coordinate. And it appears in the 4th-PN order correction in GW waveform, which requires high-precision detection of GW signal especially at lower frequencies (Seto et al. 2001;Nishizawa 2012). Fortunately, as one of its major objectives, DECIGO is competent for direct measurement of the universe acceleration.
In the future, DECIGO will detect a large number of neutron-star 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.
Summarily, the GW observations provide a powerful and novel method to estimate the spatial curvature in different cosmological-model-independent ways. 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.