The constraint ability of Hubble parameter by gravitational wave standard sirens on cosmological parameters

In this paper, we present the application of a new method measuring Hubble parameter $H(z)$ by using the anisotropy of luminosity distance($d_{L}$) of the gravitational wave(GW) standard sirens of neutron star(NS) binary system. The method has never been put into practice so far due to the lack of the ability of detecting GW. However, LIGO's success in detecting GW of black hole(BH) binary system merger announced the potential possibility of this new method. We apply this method to several GW detecting projects, including Advanced LIGO(aLIGO), Einstein Telescope(ET) and DECIGO, and evaluate its constraint ability on cosmological parameters of $H(z)$. It turns out that the $H(z)$ by aLIGO and ET is of bad accuracy, while the $H(z)$ by DECIGO shows a good one. We simulate $H(z)$ data at every 0.1 redshift span using the error information of $H(z)$ by DECIGO, and put the mock data into the forecasting of cosmological parameters. Compared with the previous data and method, we get an obviously tighter constraint on cosmological parameters by mock data, and a concomitantly higher value of Figure of Merit(FoM, the reciprocal of the area enclosed by the $2\sigma$ confidence region). For a 3-year-observation by standard sirens of DECIGO, the FoM value is as high as 170.82. If a 10-year-observation is launched, the FoM could reach 569.42. For comparison, the FoM of 38 actual observed $H(z)$ data(OHD) is 9.3. We also investigate the undulant universe, which shows a comparable improvement on the constraint of cosmological parameters. These improvement indicates that the new method has great potential in further cosmological constraints.


Introduction
In the twenty-first century, we witnessed the bloom of precision cosmology. Precision cosmology even ranked second on a list named "Insights of the decade" from Science magazine in 2010 [1]. The key of accurate cosmology is to accurately constrain cosmological parameters and their state equations, which can lead us to a better understanding of the evolution of our universe. Four main observations have been developed to constrain cosmological parameters so far [2] : Supernova(SN, [3]), Baryon Acoustic Oscillation(BAO, [4]), Galaxy Cluster(CL, [5]), Weak Lensing(WL, [6]). Actually, a relatively new tool, observational Hubble parameter data(OHD), is becoming increasingly popular these years because of its effective constraint on cosmological parameters. H(z)'s high efficiency lies on the fact that it is the only observation that can directly represent the expanding history of our universe. Compared with the Luminosity distance(d L ) of SN, H(z) contains no integral terms and directly connects with cosmological paa Corresponding author. Email: tjzhang@bnu.edu.cn rameters, which makes it powerful in constraining cosmological parameters, because the integral term can conceal many details and hide important information. As Ma and Zhang [7] reported, H(z) constrains cosmological parameters much tighter than the same-number SN does. To achieve the same constraint effect of SN subset ConstitutionT, ones need only 64 H(z) data sets under gaussian prior on H 0 , H 0 + σ H 0 .
There are various ways to detect H(z), which can be generally classified into three types: 1, differential age method [8]; 2, radial BAO method [9]; 3, standard sirens method [10]. The first two techniques have been employed in the past measurement of H(z), but the number of observed hubble parameter data(OHD) are still insufficient. We get only 38 OHD sets so far, whose accuracies are far from desirable. Now with the development of gravitational wave(GW) detecting technology, it is time to look forwards to the third method: GW standard sirens. GW standard sirens was first proposed and discussed by Schutz [10] . Schutz presented his idea that one can determine the Hubble constant through the observation of GW emitted by decaying orbit of neutron star(NS) binary system. In 2015, although the second generation GW detector Advanced LIGO(aLIGO) operated even not at its design sensitivity, it still detected the first GW signal at its first run [11]. According to theoretical understanding, GW formula of compact binary system encodes the information of d L , providing an access to the direct measurement of d L , the crucial parameter we use in this paper. We sense the possibility and potential from detecting gravitational wave.
The detection of GW not only conforms to the general relativity, but also let us see the hope of standard sirens [12][13][14]. Toshiya Namikawa et al. [15] studied GW standard sirens as a cosmological probe without redshift. But it is hard to get the corresponding redshift information without electromagnetic counterpart. The absence of redshift impedes some further research. If one were given the H(z) and its corresponding redshift, the scope of research would be wildly broaden. In 2006, a new way to narrow the relative error of H(z) by the dipole of d L has been proposed [16]. The relative error is measured by the dipole of SN. But the problem is that the new method needs plenty of SNs if we want to get a relatively accurate H(z), which can not be met in reality. Although this method also has problem in measuring high-z H(z), it is an instructive idea. Atsushi Nishizawa and Atsushi Tamga et al. [17] gave us an alternative by pointing it out that we can get information of d L through the gravitational wave function of NS binary system, instead of SN . We follow his idea and choose NS binary system as our research subject in this paper, because a rough estimate would tell us that the number of observed NS binary system turns out much bigger than that of SN. NS binary system can help us dramatically narrow down the statistical error.
Pozzo [18] proposed a general Bayesian theoretical framework for cosmological inference, which can conveniently include the prior information about the GW source. This framework defines the likelihood based on the difference between the strain of each detector and the GW template, and the posterior probability distribution for the cosmological parameters is calculated through the quasi-likelihood obtained by marginalizing over the GW signal intrinsic parameters. Applying the framework the author constrained the Hubble constant H 0 to an accuracy of 4 − 5% at 95% confidence. Nearly all subsequent work of using GW sources for cosmological inference is based on this Bayesian framework. The same framework was adopted by Taylor et al. [19], but the likelihood was defined on the assumption that the number count of GW events detected by a detector is a Poisson distributed random variate. They measured the Hubble constant using GW signals of NS binaries by narrowing the distribution of masses of the underlying NS population. That is, H 0 was determined to ±10% using ∼ 100 observations. By assuming that the masses of NS binaries can be modeled by a Gaussian distribution and that both masses of the double NS systems are equal, the authors found their chirp masses are approximately normally distributed and got the corresponding mean and standard deviation. Then, using the same method, they explored the prospects for constraining cosmology using GW observations of neutron star binaries by the proposed Einstein Telescope (ET), a third-generation ground-based interferometer. This time they fixed H 0 , Ω m,0 and Ω Λ ,0 and constrained the dark-energy equation of state (EOS) parameters [20]. With a 10 5 -event catalog, they constrained the dark-energy EOS parameters to an accuracy similar to forecasted constraints from future CMB+BAO+SNIa measurements. Chen et al. [21] investigated the measurement of Hubble constant at various cases: with and without electromagnetic counterpart, binary NS mergers and binary black hole mergers. They showed that that LIGO and Virgo can be expected to constrain the Hubble constant to a precision of 2% within 5 years and 1% within a decade. Vitale and Chen [22] dealt with neutron star black hole mergers and focused on measuring the luminosity distance to a source. They concluded that the 1 − σ statistical uncertainty of the luminosity distance for spinning black hole neutron star binaries can be up to a factor of ∼ 10 better than for a non-spinning binary neutron star merger with the same signal-to-noise ratio. Pozzo et al. [23] investigated the accuracy of the measured cosmological parameters using information coming only from the gravitational wave observations of binary neutron star systems by the Einstein Telescope. They used Fisher matrix method to extract redshift information of a source given that information about the equation of state of the source is available [24]. They found by direct simulation of 10 3 detections of binary neutron stars, H 0 , Ω m , Ω Λ , w 0 and w 1 can be measured at the 95% level with an accuracy of ∼ 8%, 65%, 39%, 80% and 90%, respectively. Different to the previous studies that focussed on constraining the parameters of specific cosmological models, our work emphasises a model independent measurement of H(z). A model free approach will generally produce a weaker constraint on any particular model than the model-specific analysis, but it has more flexibility if the true model deviates from the model assumed.
For the current observational status of GW, several frequency windows of its are targeted by different detectors. The second generation detector are mainly aimed at frequency window 10 ∼ 1000Hz, such as aLIGO and VIRGO. The next generation detector plan to reach lower frequency band. The project DECIGO was designed most sensitive at 0.1 ∼ 10Hz, while the Einstein Telescope(ET) may also reach the frequency ∼ 1Hz. The space-based eLISA can even detect GW of 10 −4 ∼ 10 −1 Hz. In this paper, we make use of GW sirens to measure H(z) by estimating the error of d L , a little different from the method proposed by Schutz [10]. Because NS binary system are used as the source of GW in this paper, the GW signal frequency of whom mostly ranges in 10 ∼ 1000Hz, we ignore the projects whose optimal sensitivity are far away from 10 ∼ 1000Hz, such as eLISA, and choose the ones whose optimal sensitivity locate around 10 ∼ 1000Hz. Finally, aLIGO, Einstein Telescope(ET) and DECIGO are chosen as our research subject.
Most importantly, the Advanced LIGO and Advanced Virgo gravitational-wave detectors made their first observation of a binary neutron star inspiral, and detected the signal of GW170817 with a combined signal-to-noise ratio of 32.4 [13,14]. In addition it provides the first direct evidence of a link between binary neutron star mergers and short γ-ray bursts. The combined analyses of the gravitational-wave data and electromagnetic emissions are providing new insights into independent tests of cosmological models, so GW170817 marks the beginning of a new era of cosmology. Using the data of GW170817, they performed the gravitational-wave standard siren measurement of the Hubble constant [14] to be 70 +12 −8 kms −1 Mpc −1 . Different from their works, in this paper, we focus mainly on two aspect: 1, How will it work out if we apply the new method to some projects? 2, how about the quality of the H(z) by this method, or to what degree could we constrain cosmological parameters? This paper is organized as follows. In sec 2, we sketch the idea of GW standard sirens method, and apply it to some GW detecting projects. In sec 3, we simulate the H(z) data, and analyze the constraint ability of the mock data for Λ CDM and the undulant universe. In sec 4, we discuss the result and talk a little about the corresponding redshift. All through this paper, we adopt the natural unit, c = G = 1.

Dipole of luminosity distance
If the universe is completely homogeneous and isotropic on large scale, and the observer is relatively rest with the cosmic microwave background(CMB), the luminosity distance, d L , would be just the same form and has the same expression as in standard cosmology. But in fact, there are perturbations around ideal condition leading into the appearance of correction term of d L [25]. Therefore d L can be written as follow: where d L represents the traditional meaning of luminosity distance in unperturbed Friedmann universe, also the average of d L on all direction, and d (1) L means the dipole of d L . The contribution to higher order terms coming from the weak gravitational lensing effect is so small compared with dipole that we ignore them here [26]. The dipole is dominated by the peculiar velocity of observers. If you want to check it or feel intrigued by the theory, you can look up the reference for the details [16].
Here is the final result: where |v 0 |, z, H(z), ∆ H(z) respectively denote the projection of observer peculiar velocity on the direction of sight, the redshift of the observed celestial body, the expanding rate at the redshift z, the absolute error of H(z), and ∆ d L respectively. From the equations above,given the value of d L , which can be solved by analyzing observed GW function in following subsection. Also, the mean error ∆ H(z) reduces to ∆ H(z)/ √ N if we observe N independent sources at the given redshift. Thus, we can improve the accuracy by the observation of a large number of sources.

GW standard sirens
One can use SN to illustrate the method of reducing the error of H(z). But due to the number and distribution of SN, it doesn't L at different redshift for v 0 = 369km/s [27]. As shown in the picture, the ratio goes very large , even bloom up, at low redshift. That is caused by the fact that the ratio approximate to (1 + z)|v 0 |/z at the limit of z = 0. work well, especially at high-z region. Considering the advantage of the larger number of observed sources, which can dramatically narrow down the error of H(z), we choose NS binary system as an alternative of SN. There is another problem for black hole binary system: black hole seldom radiates electromagnetic wave, rendering it hard to measure its corresponding redshift up to now. This is an important factor to choose NS binary system.
In GW experiments, one can extract the property of the source and cosmological information by comparing detected waveform with theoretical template. That is what LIGO team did when the first detected the GW of two back holes merged [11]. The typical Fourier transform of GW waveform can be expressed by which is based on the average over sky location. Here A = ( √ 6π 2/3 ) −1 is a constant geometrically averaged over the inclination angle of a binary system. d L (z) is the luminosity distance at redshift z, and we can set it as d (0) L because we need to observe plenty of source at the given redshift. M z = (1 + z)η 3/5 M t with the definition of total mass M t = m 1 + m 2 and symmetric mass ratio η = m 1 m 2 /M 2 t . The last unknown function Ψ ( f ) is a little intricate. It is the frequency-dependent phase caused by orbital evolution. Usually we deal with it by post-Newtanion(PN) approximation, an approximation to general relativity in the weak-field, slow-motion regime [28]. Its concrete expression will not affect the final result, because this term will be eliminated when we do the following calculation.
Here we just need to know that it is a function of the coalescence time t c , the phase when emitted φ c , M z , f , η.
There are five unknown parameters, namely: M z , η, t c , φ c , d L , where d L is the only parameter that has nothing do with the own property of binary system. For the convenience of calculating, we just take account of equal mass NS binary system with 1.4M , and set t c = 0, φ c = 0. Then we have M z = 1.22(1 + z)M , η = 0.25. Though GW may tell us some information about the redshift [24,29], we have no data about the redshift and we need a general method to get the reshift information. We should still resort to electromagnetic observation to find out corresponding redshift. Cutler and Holz [30] demonstrated its technological viability. Fig. 2: The noise power spectrum. Green curve represents P 1 ( f ) for DECIGO, blue curve represents P 2 ( f ) for ET, red curve represents P 3 ( f ) for aLIGO respectively.
We use fisher matrix to estimate error. Fisher matrix has its limit: the Cramer-Rao bound, and it's break down at low SNR. The error estimate for d L is based on Fisher matrix that is given by where ∂ a means derivative with respect to parameter θ a . For DECIGO, there exist eight interferometric signals, Γ ab should multiplied by 8. We set values to parameters expect for d L , so the only parameter in Γ ab is d L . P( f ) is the noise power spectrum, and the P( f ) for DECIGO, ET and aLIGO are shown in Fig. 2. Here we give the expression of each detector's noise curve, P 1 ( f ), P 2 ( f ), P 3 ( f ) for DECIGO, ET and aLIGO respectively.
DECIGO: DECIGO(Deci-Hertz Interferometer Gravitational Wave Observatory) is a planed space-based GW observation aimed at 0.1 ∼ 10Hz frequency window. Its configuration is still to be decided. Here we adopt the following parameters in its configuration: the arm length 1000km, the output laser power 10W with wavelength λ = 532nm, the mirror diameter 1m with its mass 100Kg, and the finesse of FP cavity 10. Thus its noise curve is [31] ET: ET(Einstein Telescope) is a third generation GW detector, whose design is not finished. Here we just consider the simplest case with 10km arms. We adopt the fitting expression given by Keppel and Ajith [32] aLIGO: aLIGO is an available second generation detector whose optimal sensitivity band match with the frequency window of GW from NS binary system . The first run of aLIGO did not reach its design sensitivity. Here we use the noise curve fitted by [33]. It is not an accurate expression, but an approximation of the original curve is given by [34] In the expression of Γ ab , the lower cutoff of frequency, f min , is a function of observation time T obs , f min = 0.233(1M /M z ) 5/8 (1yr/T obs ) 3/8 In the case of our paper, for a given T obs , f min changes little with M z , which is always in the high strain noise region. It makes no big difference to the result of the integral. A reasonable setting of the value of f min will work. But for prudence, we just take the original expression of f min when calculating the integral. And the higher cutoff, f max , is the inner-most stable circular orbit frequency, whose typical value is of kHz order [35]. More specific, f max 2000Hz in our case. In the calculation, f max can be set by the property of the integrand. The value of integrand sharply drops with f getting larger, so its contribution to Γ ab can be neglected. For the reason of integrand property, we set the f max of DECIGO, ET, aLIGO respectively to be 100Hz, 2000Hz and 2000Hz. Then the one-sigma instrument error is If we launch an observation for a given source, the one-sigma error estimate σ instr of d L arises from instrumental noise. For a given device, no matter it is DECIGO, ET or aLIGO, the accuracy of d L is the same even for different observation time. It makes no difference for the σ instr no matter how long the observation continues, which is mainly because that the error is For any observation longer than that time the precision is the same since you do not observe the source any more. The σ instr of three devices are showed in Fig. 3. As we can see, the accuracy is far from desirable. The H(z) error would increase if we include other errors . We need to take measure to narrow down the error. This is what we do in next subsection

H(z) error
In last subsection, we already calculate the device error σ instr for a given NS binary system under a given observation device.
Besides the device error-the dominating error, there are two main errors, namely the lensing error and the peculiar-velocity error. The lensing error is due to the lens effect. Here we take a recent fitting by [36], And the peculiar-velocity error is a kind of Doppler effect of the movement of the celestial body, essentially. The peculiarvelocity error can be described as ([37]) where σ v,gal = 300km/s is the approximation of the 1-dimensional velocity dispersion of the galaxy. Then we get the expression of relative error of d L (z): Before we do the calculation to get the relative error of H(z), there is one more step we can do for a better accuracy. The mean error will statistically abate if we have many independent sources. Reducing the error of H(z) by observing many NS binary system at the same redshift may be feasible. The problem is to what degree can we reduce the error? First we need to figure out the number distributionṅ(z) of NS binary system at different redshift. The distribution of NS binary system can be described and estimated. According to Cutler and Harms [38], the fitting of NS-NS merger rate can be given by: where s(z) is estimated from star formation history inferred from UV luminosity, andṅ 0 represents the merger rate at present time. Then ∆ N, the number of NS-NS merger at redshift bin ∆ z, is expressed by: ∆ N(z) = T obs z+∆ z/2 Recent work doesn't provide solid evidence of the exact value ofṅ 0 . he latestṅ 0 range inferred by the observation of GW170817 is 0.32 − 4.7 × 10 −6 Mpc −3 yr −1 [13]. Also not every merger event would be detected. Here we encounter an conundrum. Considering that we are aimed at evaluating the method, not launching an actual observation here, we decide to, a bit arbitrarily, setṅ 0 equal to 1.0 × 10 −6 Mpc −3 yr −1 , and assume all the merger events could be detected, and the redshift width ∆ z = 0.1. Thus we get the estimate of 10-year observed number of NS binary system merger at different redshift, which is shown in Fig. 4.
The total number of SN is just of hundred-magnitude by now, while the observed number of NS-NS merger event will be much larger than that of SN, showing a tremendous potential in reducing the mean error of H(z). And from above equation, the number of NS-NS merger at fixed redshift increases with T  Combining with the information we get above, we can calculate the H(z) error for a specific device under a given observation time now. The relative error of H(z) by DECIGO, ET and aLIGO is shown in Fig. 5, Fig. 6, and Fig. 7, for 1-year, 3-year, 10-year observation respectively. The relative error by aLIGO is a total disaster, which basically has little application value in constraining cosmological parameters. The error by ET is a little better, especially at low redshift region, because ET is more sensitive than aLIGO. Thus DECIGO plays best in this method. When redshift reach 3, due to the decreasing of the number of observed NS-NS merger event with redshift, the relative error of H(z) is magnified, but still quite small. And the elongation of T 1/2 obs shows a great ability in narrowing down the error. We stress here that σ lens (z) contributes a lot to the distance error ∆ d L . Of course, various techniques have been developed to reduce σ lens (z). Stefan Hilbert [39] suggested that deep shear survey can narrow the lens error. Hirata [36] found that ones can improve the distance determination typically by a Fig. 7: Same as Fig. 5, but for aLIGO. For a given observation time, the dash line and the solid line overlaps, because the lens error is relatively small compared with the instrument error of aLIGO factor of 2-3 by exploiting the non-Gaussian nature of the lens magnification distribution. C. Shapiro [40] used the procedure 'delensing', to estimate the magnification and thereby remove it by a weak lensing map. It may be too optimistic to remove all the lens error. But we can rely on it that we could reduce the lens error to a very low level in the near future. Therefore, for simplicity, we will ignore σ lens (z) in the following sections.

Simulating data
The new method for measuring H(z) has been presented and its error analysis has been done above. The problem is how H(z) can accurately observed by this way to constrain cosmological parameters? So far, we did not obtain actual OHD by this way. But it doesn't necessarily mean we can do nothing about it. A reasonable and rational simulation can help us forecast and evaluate. Since H(z) by aLIGO has a bad accuracy, we carry on no simulation and forecast for aLIGO here. H(z) by ET can do some simulation and forecast. The problem is that the effect is a little bad, even worse than 38 OHD sets. We do not plan to show it here, too. Thus, DECIGO is the only device we discuss in following sections. Now that we have the error information of H(z), we can simulate the OHD. We follow the method Yuan Shuo [41] to generate mock data for Λ CDM:  Fig. 8. We get 38 OHD sets up to now. The datas were obtained by different ways from different groups [8,9,[42][43][44][45][46][47][48][49][50][51]. Fig. 9 shows H ΛCDM at every redshift and the 38 OHD sets so far. As we can see, the OHD value goes up and down around the H ΛCDM at the same redshift, which justifies the validity of our simulation. Compared with our mock data, the actual OHD sets' accuracy is obviously worse.

Forecasting
Now that we have got the 3-year-observation mock data, we can use them to forecast. Before that, we need a criteria to evaluate the constraining ability of the dataset-Figure of Merit(FoM).
We can define FoM in different ways, as long as its value can reflect how tightly or loosely the data constrains parameters.
Here for the convenience of our analysis, we adopt the definition by Albrecht [52], the reciprocal of the area enclosed by the 2σ confidence region contour, coinciding with a specially appointed confidence region under gaussian distribution.
We choose the Λ CDM as our prior model. In such a standard Λ CDM universe with a curvature term Ω k = 1 − Ω m − Ω Λ , the Hubble parameter is given by The determination of H 0 has been carried on in different H 0 tension projects. For 7-year WMAP observation, H 0 = 73 ± 3kms −1 Mpc −1 [53]. In this paper, we take the most recent value H 0 = 74.2 ± 3.6kms −1 Mpc −1 [54]. The best value of Ω m , Ω Λ we adopt is 0.27, 0.73 respectively, due to the coherence that they are consistent with the observations and the fact that we use these value to generate our simulation data. All the three parameters are assumed under gaussian distribution. By Bayes' theorem, the posterior probability density function of parameters is where is the likelihood and P(H 0 ) is the prior probability density function of H 0 . And the expression of is given by(assuming no covariance between parameters) where σ i is the uncertainty of the data H i . And the P(H 0 ) is Gaussian prior, given by Then the integral can be worked out for the given Gaussian prior P(H 0 ). There is a point in the parameter space maximizing the probability density, P max . Because of what we have described in last paragraph, such a point in this forecasting is {0.27, 0.73, 74.2}. The formula means the contour of a given confidence region, which corresponds to the value of ∆ χ 2 . We have three parameters here, Ω m , Ω Λ , H 0 . ∆ χ 2 is statistically set to 2.3, 6.17, 11.8 respectively for 1σ , 2σ , 3σ confidence region. For a direct comparing and understanding, here we choose 2σ confidence region, namely ∆ χ 2 = 6.17, when calculating FoM. To order to calculate the confidence region and FoM, we take the Fisher Matrix forecast technique [55], where the value of matrix elements is taken at the most-likely value of parameters. Let's denote the marginalized Fisher matrix byF, then the contour in subspace is given by where ∆ θ is the deviation from the beat value of the parameters. When calculating FoM, we take ∆ χ 2 as 6.17. The enclosed area is π/ det(F/∆ χ 2 ). So FoM, the reciprocal of the area, is The contour is shown in Fig. 10. As we can see, the contour is an ellipse, which is consistent with the equation ofF. For a more direct and concrete comparison, we perform constraint for the 38 OHD sets. Their constraint on Ω m and Ω Λ is shown in Fig. 11. Apparently, the constraint of the mock data on parameters is much tighter, compared with that of available OHD sets, which has an significant improvement on precision cosmology. The simulation and forecast of 10-year-observation is just carried out in the same way. As Fig. 12 shows, its constraint on cosmological parameters is even much tighter, implying a consequent higher FoM value.
For the FoM value of 3-year-observation mock data, it is about 170.82, while that of 38 OHD sets is just about 9.3. It is a remarkable improvement. For 10-year-observation mock data, the FoM has a farther improvement, reaching 569.42. We have enough reason to look forward to the excellent application of H(z) data by this method.

Nonstandard model
The Λ CDM universe do match the observation quite well. But it doesn't answer the question that why matter and vacuum energy should be of the same order of magnitude at this moment. Here we consider another model which can give us an answer to this problem by alternating periods of acceleration and deceleration. In undulant universe, the equation of state of the vacuum energy is an oscillatory function of state of the scale of the universe, w(a) = −cos(Ina). It meets the fact that w(a = 1) = −1 in the current universe. Then the Hubble parameter is given by: where Ω Λ +Ω m +Ω k = 1. The simulation and forecasting carry out just the same as above. Here we consider the case of 3-year-observation. The corresponding FoM is 153.07. And the constraint is displayed in Fig.13

Discussion and Conclusions
In this paper, we mainly evaluate the quality of H(z) data by GW standard sirens method of several GW detection plans, whose optimal frequency locate around the frequency window of GW from typical NS binary system. We calculate the relative error of H(z) by three devices, DECIGO and ET and aLIGO. Though the sensitivity of the three devices is almost of the same order of magnitude, the H(z) error of DECIGO is quite optimistic while that of other two is far from satisfying. But it does not mean that H(z) data by this method is a dead end or of no meaning, which is justified by the forecasting of DECIGO-based H(z) data. If the sensitivity of aLIGO or ET is sightly improved, or just move the most sensitive frequency to a lower region, the error of H(z) will be comparable that by DECIGO. To get a better H(z) data, we have two ways: (1) the noise curve could be moved down, and the signal-to-noise ratio of a given sore increases, so we get more events above threshold and more "bright" events which have smaller errors; (2) the noise curve could be moved to the left, then we can see more of the inspiral which can help improve parameter estimation at fixed signal-to-noise.
Considering the absence of real H(z) data by DECIGO, we simulate H(z) and the data show an alluring constraint ability on cosmological parameters. After all, we are aimed at evaluating the viability and quality of H(z) data by GW standard sirens method, not putting the method into actual operation. We find that, under Λ CDM universe, the FoM of mock data shows a huge improvement compared with that of 38 actual OHD sets. For contrast, the FoM is 9.3 for 38 OHD sets, 170.82 for 3year-observation, 569.42 for 10-year-observation respectively. The advantage of H(z) is that it is the direct measurement of the expansion history, so H(z) can be powerful in constraining nonstandard universe. Besides the standard model, we also explore its ability when applied to undulant universe. H(z) by DECIGO still shows a excellent constraining ability and a comparably excellent result. For 3-year-observation simulation, the FoM is 153.07. The tight constraint of mock data and the FoM of the corresponding contour indicate a bright future of measuring H(z) data by this method. To extract as much physical information as possible, all the known sources of error should be quantified. Apart from the three error mentioned above, there is another kind of errors, the calibration error. In GW detection, the response function is used to convert the electronic output of a GW detector into the measured GW signal. The calibration error is produced on the experimentally measurement of the response function [56]. The calibration error in the response function degrades the ability to measure the physical properties of the GW source. Thus it is meaningful to investigate the calibration error. Lee Lindblom [56] derived the optimal calibration accuracy: the lower accuracy level would reduce the quantity and the quality of the scientific information extracted from the data, and the higher accuracy would be made irrelevant by the intrinsic noise level of the detector. And S. Vitate at al. [57] also investigated the effect introduced by calibration error based on the estimates obtained during LIGO's fifth and VIRGO's third science runs. They found that the calibration error would slightly damage the parameter estimate in GW data analysis. But the calibrationintroduced system has a better ability in locating the source, facilitating the EM counterpart detecting. Considering the damage caused by calibration error is relatively small and its hard to quantify the calibration error, we ignore it in this work. We expect future study can give us more precise estimate.
It is well worth to point out that at current stage, the detecting of electromagnetic counterpart is still a problem. Currently, we are not able to make a good evaluation and conclusion about it. Finding the EM counterpart in the GW event is crucial for GW astronomy, which can reveal the process and interaction during the merger process [58]. Mwtzger et al. [59] showed that the transient EM counterpart can possibly occur within a few seconds after the binary merger. And a lot of theoretical and experimental progress have been achieved [60,61].
Also there are some recent development in astronomical and computing technologies. During the proposal and test of a number of low-latency GW trigger-generation pipelines, the pipeline has been improved and capable of generating event triggers within minutes upon the arrival of a detectable signal [60,[62][63][64]. More and more detectors to be constructed can form a network, rendering it likely to improve the localization efficiency [65][66][67]. Some methods have been proposed to identifying GW source for a large sky error [68,69]. Considering the fact that the early detector networks error in GW localization will be of order 200deg 2 [70], such method would improve the feasibility of EM detector a lot. The devices that aim at facilitating the prompt EM detection mainly focus on high energy region and the optical region, while radio region is also a good candidate [58]. By next decade, the Large Synoptic Survey Telescope(LSST), will be in its sky survey. It will bring us great hope to find the prompt EM counterpart. Such EM detection demand multi-wavelength programs by sensitive telescope capable of covering large areas on the sky, and a strong synergy exists between LSST and radio survey in identifying the EM counterpart at both optic and radio wavelengths, and the information from both wavelengths about the physics of the postmerger will be complementary [71]. Here we stress the evaluation of H(z) by standard sirens, not the exact technical details. Another problem is that the detecting range for NS binary system is just 300Mpc now [72]. This range is much smaller than what we assumed above. We explore how the FoM changes with the variation of detecting range, which is shown in Fig. 14. The FoM can be comparable with that of 38 OHD sets at z = 0.7. For 10-year-observation, this critical value would be smaller.In other words, if we launch a 10-year-observation, even if the detecting range is just z=0.7, we can do much better than 38 OHD sets. The reason why the limited data can produce such good effect lies on the fact that the GW standard sirens can measure low-z H(z) with excellent accuracy. This demonstrates that even if we could not detect the high-z H(z) data by GW standard sirens method, the low-z data can still be valuable and powerful. In the further, if we want to measure high-z H(z) by this way, some improvement, probably a lower strain noise, is necessary. But it is undoubtable that the H(z) by this method is of great power and potential.

acknowledgements
We sincerely thank the referee for his/her useful responses, which help us greatly improve our manuscript. This work was supported by National Key R&D Program of China (2017YFA0402600), the National Science Foundation of China (Grants No. 11573006, 11929301, 61802428) and Shandong Provincial Natural Science Foundation of China(Grant No. ZR2019MA059).