How will our knowledge of short gamma-ray bursts affect the distance measurement of binary neutron stars?

GWs from BNS associated with SGRBs have drawn considerable attention due to their prospect in cosmology. For such events, the sky locations of sources can be pinpointed with techniques such as identifying the host galaxies. However, the cosmological applications of these events still suffer from the problem of degeneracy between luminosity distance and inclination angle. To address this issue, a technique was proposed in previous study, i.e., using the property of SGRBs. Based on the observations, we assume that the cosine of inclination follows a Gaussian distribution, which may act as a prior in the Bayes analysis to break the degeneracy. This paper investigates the effects of different Gaussian priors and detector configurations on distance measurement and cosmological research. We first derive a simplified Fisher information matrix for demonstration, and then conduct quantitative analyses via simulation. By varying the number of third-generation detectors and the scale of prior, we generate four catalogs of 1000 events. It is shown that, in the same detecting period, a network of detectors can recognize more and farther events than a single detector. Besides, adopting tighter prior and employing multiple detectors both decrease the error of luminosity distance. Also considered is the performance of a widely adopted formula in the error budget, which turns out to be a conservative choice in each case. As for cosmological applications, for LCDM model, 500, 200, 600, and 300 events are required for the four configurations to achieve 1% H_0 accuracy. With all 1000 events in each catalog, H_0 and Omega_m can be constrained to 0.66%, 0.37%, 0.76%, 0.49%, and 0.010, 0.006, 0.013, 0.010, respectively. The results of the Gaussian process also show that the GW standard siren can serve as a probe of high-redshift universe.


Introduction
Ever since the first observation of gravitational wave (GW) emitted by binary neutron star (BNS), known as GW170817 [1], and the coincident short gamma-ray burst (SGRB) event GRB170817A [2], the method of using gravitational wave standard siren (GWSS) [3] together with electromagnetic (EM) counterpart [4] to study cosmology has come into reality. The GWs from compact binaries provide *Corresponding author (email: lxxu@dlut.edu.cn) direct measurements of the sources' luminosity distances. The redshift of a GW source can be independently obtained via extra information from EM counterparts [5] or with statistical methods, such as the dark siren [6][7][8]. Among various EM counterparts, SGRB, which is usually associated with BNS, appears especially useful. By identifying the host galaxy of SGRB, it is possible to determine the redshift and sky position of a GW source at the same time. Moreover, the collimation property of SGRB can also give a clue to the source's orientation. Although GW sources of other types The cosmological applications of GW signals rely on determining source parameters, especially the luminosity distance. Intrinsic parameters, such as the component masses, can be measured with exquisite accuracy from the phase of GW [20]. In the presence of counterparts such as SGRB, it is possible to precisely probe the sky position of the source by identifying the host galaxy. Unfortunately, since luminosity distance degenerates with the inclination angle, using GWSS as an absolute distance probe requires a mechanism to break the degeneracy [21]. Employing a network of detectors contributes to decreasing the uncertainty of distance [22]. In recent years, a customary practice is to neglect the distanceinclination correlation and account for the influence of correlation by multiplying a factor of 2 [23]. We argue that the validity of this method should be evaluated by quantitative analysis.
Several recent studies have shown that information on SGRB signals may help break the degeneracy. The authors of [24] applied a Bayesian framework to 1,000 simulated joint detections by advanced LIGO and advanced Virgo to demonstrate that combined SGRB and GW observations could improve the estimations of progenitor distance and inclination angle. Besides, investigated in [25] is the role of EM measurement in the analysis of GW data to improve the precision of luminosity distance, and hence that of the Hubble constant. This approach was tested on GW170817 and simulated events observed by the HLV and HLVJI networks. Another demonstration of the idea can be found in [26]. Based on these pioneering studies, our research is dedicated to investigating how knowledge of SGRB affect the measurement of luminosity distance, and its further influence on cosmology, with both model-dependent and model-independent methods.
In this paper, we adopt a method inspired by [21,24,25]. To solve the problem of degeneracy, we employ the propertie of SGRB in the form of a Gaussian prior on the inclination angle. The posterior distribution of parameters is then calculated using the Fisher information matrix (FIM). In the simplified regime where there is only one stationary interferometer, the FIM of luminosity distance and inclination can be analytically deduced, giving intuitive explanations to the problem of degeneracy and how it is broken by the prior. By varying the number of 3G GW detectors and the scale of prior, four catalogs of GW events are simulated to test the abovementioned pipeline of parameter estimation.
We give a detailed description of our methods in sect. 2, including the detector response of 3G GW detectors, the manipulations of FIM based on SGRB prior, and the settings of simulations. Presented in sect. 3 is the error analysis of simulated GW catalogs, followed by the application of mock data to the ΛCDM model and Gaussian process (GP). We further study the effects of altering the form of prior and star formation rate (SFR) model in sect. 4. The concluding remarks are given in sect. 5. In this paper, we adopt the natural units G = c = 1.

Gravitational wave detector response
The detector response of GWs from coalescing compact binary is briefly summarized in this section. We calculate the coefficients of vectors and tensors with respect to an imaginary geocentric frame. Given a binary system at sky position n =n (θ, φ) with the orientation of the angular momentum L, a pair of axes can be constructed as follows: wheren,L are both unit vectors, and we define the inclination angle as ι = arccos n ·L . The pair {X,Ŷ } spans the projected orbital plane orthogonal ton, andX andŶ are along the major and minor axes, respectively. By denoting the angle betweenX and the orbit's line of nodes as the polarization angle ψ, the coefficients ofX andŶ can be explicitly calculated [27]. The GW tensor in transverse traceless (TT) gauge is the sum of two polarizations where e +,× are the basis tensors of polarization, defined as Under incoming GW, the output of a ground-based Michelson interferometer (labeled I) is a timeseries: where the detector tensor of an interferometer with armsl and m is D I = l ⊗l −m ⊗m /2, and the time-dependence of D I comes from the rotation of the Earth. For simplicity, we assume that the Earth is a perfect sphere with radius R and rotational angular velocity Ω. For an interferometer at r =r(π/2 − φ I , λ I = λ I0 + Ωt) (λ I0 and φ I are the east longitude and north latitude, respectively), whosel arm is oriented at angle γ I north of east andm arm at γ I + ζ I , the unit vectorŝ l,m can be expressed aŝ where e E I and e N I are the unit vectors pointing to the east and north on the Earth's surface. Combining eq. (2) and eq. (4), the time-domain detector response reads where F +,× I = Σ i j D I,i j e +,× i j are the antenna pattern functions. A coalescing BNS system is characterized by component mass (m 1 , m 2 ), luminosity distance d L , time and phase of coalescence (t c , φ c ); thus, the mass ratio is η = m 1 m 2 /(m 1 + m 2 ) 2 , and the redshifted chirp mass is M c = (1 + z)(m 1 + m 2 )η 3/5 .
The orbital frequency varies negligibly over a single GW cycle during the inspiral section of BNS coalescence, which makes it possible to compute the Fourier transform of h I with the stationary phase approximation [17]. In this paper, we adopt the restricted post-Newtonian (PN) approximation and calculate the waveform to the 3.5 PN order [17]: where the 2π fn ·r term is included to account for the time of GW propagating from the Earth's center to the detector. Other terms are defined as follows: where d L and ι are the luminosity distance and inclination angle, respectively, and M represents the redshifted total mass, defined as M = (1 + z)(m 1 + m 2 ). Detailed expressions for the PN coefficients ψ i can be found in [28]. Notably, the time-dependence of F +,× I should be converted to frequency-dependence through F +,× . This principle also applies tor.
Three 3G detectors will be considered in the following research: ET in Europe, CE in Idaho, USA, and an assumed CE-like detector in New South Wales, Australia. The parameters {λ I0 , φ I , γ I , andζ I } are listed in Table III of [29] and Table 1 of this paper. Besides, we choose the sensitivity curves of ET and CE as ET-D and CE2-40-CBO, whose definitions can be found in [29]. We denote the network of ET and two CE-like detectors as ET + CE.

Fisher information matrix
The detector response of GW entirely depends on nine parameters: {t c , φ c , η, M c , θ, φ, ψ, ι, d L }. Alternatively, one can use GW signal to constrain these parameters, among which we are particularly interested in d L . Suppose that we have a set of observational data D and a model H with N parameters dubbed λ. According to the Bayesian theorem, the posterior of λ is where p (D|λ, H) denotes the likelihood function, and p (λ|H) is the prior. The normalization factor p (D|H) can be calculated by requiring d N λp (λ|D, H) = 1. If we choose flat priors on all parameters, then the posterior is proportional to the likelihood function.
In general, the result of full Bayesian analysis is more precise and rigorous than that of FIM, and the Bayesian method is already adopted in topics such as GWSS [30] and dark siren. While, it has been shown that the results of FIM agree with sophisticated Bayesian parameter estimation in the condition of high signal-to-noise ratio (SNR) [31,32], which, as will be seen later, is always satisfied in our simulation. The implementations of FIM in the context of GW parameter estimation can be found in recent research such as [33][34][35]. More importantly, we are especially concerned about the relationship between luminosity distance and inclination, and the simplified FIM of these two parameters can provide intuitive explanations to the problem of degeneracy and how it is broken by the prior. Therefore, in this paper, we will implement FIM as the primary approach of parameter estimation. Following [22], the joint posterior distribution of parameters can be approximated by a multidimensional Gaussian form: where is the FIM of a detector network comprising N d interferometers (N d = 3 for ET, N d = 5 for ET + CE) and ∆λ represent the deviations of λ to their "true" values. The inner product of any two functionsã( f ) andb( f ) is S h,I ( f ) being the one-sided noise power spectrum density (PSD) of the Ith detector. We set f lower = 1Hz and f upper = 2 f LSO = 1/(6 3/2 πM), sufficient to cover the inspiral section of BNS coalescence, where M = (1 + z)(m 1 + m 2 ) is the redshifted total mass (see [36][37][38] for the PSDs of ET and CE). We will find it convenient to give the definition of SNR: According to eq. (7), Once we have Γ i j , it is straightforward to evaluate the accuracies of parameters, characterized by the covariant matrix Σ = Γ −1 . Thus, the root-mean-square (RMS) error of λ i is σ ∆λ i = √ Σ ii , and the covariance between two parameters λ i , λ j is c i j = Σ i j / Σ ii Σ j j .
In the context of BNS-SGRB event, the sky position (θ, φ) of the GW source can be pinpointed by techniques such as identifying the host galaxy. Moreover, the mock data challenge of the 3G GW detector [39] indicates that the mass parameters can be constrained to great accuracy. Since our interest is on the determination of luminosity distance, parameters such as t c and φ c , which only appear in the expression of phase, only have negligible impact. These considerations leave only (ψ, ι, d L ) to be determined. Further investigations on the parameter space [21,22,40] have shown that the correlation between d L and ι is overwhelmingly intensive; thus, we focus only on the covariant matrix of d L and ι and assume that other parameters are perfectly constrained.
It then comes down to calculating the derivatives ofh I ( f ). Below are some relevant expressions, and the derivatives with respect to other parameters can be found in [41,42]: where we have replaced the parameter ι by v = cos ι for convenience, and the real and imaginary parts of ∂ lnh I /∂v are abbreviated as v 1,I and v 2,I , respectively. Substituting eq. (18) and eq. (19) into eq. (14), we have Generally, v 1,I and v 2,I are frequency-dependent due to the rotation of Earth, and fixing the Earth's orientation would increase resulting errors. However, the form of FIM would become quite simple under this assumption: where the elements are arranged in the order (d L , v). Further, in the case where there is only one interferometer (N d = 1), dropping the suffix I, the inverse of FIM is as follows: Next, we will illustrate our technique of handling FIM under these simplifications and perform qualitative analyses. Whereas the full expression eq. (20) will be used in the simulation (sect. 2.4). The SGRBs are believed to be strongly beamed phenomena [43][44][45], which was confirmed by recent investigations on GRB 170817A. From theoretical expectation, it is usually assumed that for SGRB induced by BNS coalescence, the viewing angle of SGRB is identical to the inclination angle ι. [46] proposed that ι follows Gaussian jet profile with mean value 0 and standard deviation σ ι = 0.057 +0.025 −0.023 rad, while [47] constrained σ ι = 4.7 +1.1 −1.1 deg. In other previous research [21,40], it was assumed that SGRB is confined within 25 • . As a result, the behavior of Σ near v = cos ι → 1 (face-on) is the most relevant.
Since face-on and face-off are equivalent, we focus only on positive values of v. It is easy to figure that when v → 1, v 1 → 1 and v 2 → 0; thus, the errors of d L and v (dubbed σ ∆d L and σ ∆v ) diverge, and the correlation c d L v = v 1 / v 2 1 + v 2 2 approaches 1, meaning that d L and v are intensively degenerate. Moreover, in the limit of v = 1, Γ is singular so that Σ cannot be well-defined, despite the fact that SNR actually increases for larger v.
To address this issue, several methods have been proposed. In recent years, a customary practice in the distance measurement of GW, as depicted by [23], is neglecting the correlation between d L and v, calculating the error of d L at v = 1 as d L /ρ, and accounting for the influence of varying v with a factor of 2 (i.e., σ ∆d L = 2d L /ρ, dubbed σ ∆d L ,0 hereafter) [38,[48][49][50]. We argue that this is different from considering the correlation at the first place; thus, the validity of using σ ∆d L ,0 should be assessed.

Short gamma-ray burst prior
In this section, we attempt to solve the problem of degeneracy by considering a nonuniform prior on v, assuming that we have knowledge about the properties of SGRB before detection. As stated in sect. 2.2, SGRB should be preferentially beamed along the orbital momentum axis of BNS. Following [21,40], the prior on v takes a Gaussian form with standard deviation σ v : This is a pessimistic scenario where the inclination angle is determined with an accuracy of, e.g., ∼5 deg at 1σ [46,47]. Moreover, σ v is a global parameter so that a single value can be used to describe a whole population of astrophysical events that are uncorrelated. Combining eq. (12), eq. (13), and eq. (23), the posterior distribution of parameters is as follows: where ∆λ = (∆d L , ∆v).
In general, the distribution given by eq. (24) is non-Gaussian due to the "true" value of v is usually not exactly 1. However, we do not expect the errors to deviate from the face-on case considerably. As a schematic demonstration, we first analyze the situation where the "true" value of v is 1. In that case, p(∆λ) returns to the two-dimensional Gaussian form. Therefore, we can still describe the parameter accuracies in the framework of FIM. For the simplified regime described in sect. 2.2, we have Interestingly, σ ∆d L /d L = 1/ρ 2 + σ 2 v gives an estimation comparable to σ ∆d L ,0 = 2d L /ρ, and the catastrophic error near v = 1 is avoided. Besides, d L is better constrained when ρ is larger or σ v is smaller, consistent with intuition. Also notable is that the correlation In the simplified regime, we further consider the situation where v 1 and investigate the effect of varying v. To this end, a bunch of events detected by an interferometer with ET-D sensitivity is simulated, with the mass parameters of GW170817 (M c = 1.188M sun , m 2 /m 1 = 0.85), located at redshifts 0.5, 1, 2, and 5. For reasons given below, two values of σ v (0.05 and 0.003) are adopted. When v 1, eq. (25) cease to be valid, and the FIM depends on F +/× , hence the angles θ, φ, ψ; thus, we average σ ∆d L over these parameters. Shown in Figure 1 are the relative errors of d L as functions of v, with v ranging from 0.75 to 1, equivalent to the 5σ range for σ v = 0.05. The solid and dashed lines correspond to the results of σ v = 0.05 and σ v = 0.003, respectively. In each case, the error of d L decreases as v approaches 1. Moreover, with tighter prior on v, one can estimate d L more accurately, and the influence of σ v becomes insignificant for distant events whose SNRs are relatively lower.
Notably, eq. (25) merely serves to provide a qualitative illustration, and we will use the full expression eq. (24) in simulation. The error of luminosity distance is calculated by integrating over the entire distribution, following the standard definition in statistics.
In this paper, two values of σ v will be considered: σ v = 0.05, consistent with [21,40] and σ v = 0.003, equivalent to the 1σ range of [47,48]. The first value seems rather large compared with the observation of GRB 170817A [47], yet it can act as an upper bound for future observations. Although future SGRB observations may give σ v other than these two values, our research give a clue on the changes of results when σ v varies.
The overall performance of our method depends on the distribution of source parameters. Details are presented in sect. 2.4.

Simulation
The fiducial cosmological model is flat ΛCDM with parameters H 0 = 67.8km s −1 Mpc −1 and Ω m = 0.308; thus, we have The redshifts of sources are drawn from the distribution where the second equality comes from the definition of comoving volume, and H(z) is the Hubble parameter at redshift z. The merger rate of BNS R merge (z) is calculated in the source frame, which is the convolution of neutron star (NS) formation rate R f and a time-delay distribution P d (t d ) [48,51]: where t d is the time between the formation and merger. By neglecting the lifetime and mass variation of the progenitor star [52], R f follows the SFR with parameters of the "Fiducial + PopIII" model [53]. Further, we assume P d (t d ) ∝ t −1 d for t min = 20Myr and t max = t H , t H is the Hubble time.
As for other parameters, the sky locationn =n(θ, φ), polarization ψ, time and phase of coalescence (t c , φ c ) are sampled from uniform distributions, and v is drawn from eq. (23). The component masses of BNS follow the Gaussian distribution N(1.33, 0.09)M sun , according to [54]. Besides, we classify a GW event as detectable if the total SNR exceeds the threshold of ρ threshold = 12 for ET + CE or 8 for ET. We set a lower threshold for ET to not limit its observation depth, making the population detected by ET and ET + CE comparable. These thresholds are widely adopted in the literature.
Notably, current SGRB observations accumulate at only low redshifts, which is partly because the apparent luminosity of SGRB decreases with distance, and an event has to exceed the flux limit of the detector to be detected. However, from the theoretical perspective, coalescing BNS can be accompanied by SGRB provided proper conditions are satisfied. Therefore, we optimistically assume that SGRB counterparts at redshifts higher than the range of current observation can be detected in the future. We believe that with a reliable prediction of the SGRB detection rate, the analysis of our study would become more practical, while this is beyond the scope of this paper.
We do not intend to predict the exact detection rate of BNS-SGRB events, given that we only have very few observations to date. Contrariwise, we will generate a catalog comprising 1,000 events (for each configuration), and for every single event, the detecting criteria of both GW and SGRB are satisfied. This requires that more than 1,000 sources should be simulated, and then the rate of acceptance is defined as 1,000 divided by the number of all (qualified or unqualified) sources. For current SGRB observation by Swift, the number of SGRBs with redshift measurement is less than 10 per year. However, the authors of [48] simulated GW detection with an ET + CE + CE network, assuming that a THESEUS-type [55-57] satellite will be used for coincidence searches, and they obtained 907 GW-SGRB joint detections in 10 years of data collection if the mass distribution of BNS is Gaussian (as is assumed in our study). Therefore, we consider it reasonable to expect 1,000 joint events with the help of future THESEUS telescope. Next, we will study how many events are sufficient for the target accuracy required by cosmological research. Four configurations will be considered in this paper:

Simulated catalogs
In the former section, four catalogs of 1,000 GW events were simulated. Shown in Figure 2 are the redshift distributions of sources, where the curves labeled "theory" are plotted according to eq. (27). It is worth noting that the histograms of A and B fit the theoretical curve well. Indeed, up to 78.9%(A) / 82.4%(B) of the entire samples pass the SNR threshold of ET + CE.
On the other hand, data in catalogs C and D tend to distribute at lower redshifts, since SNR is inversely proportional to d L , and the SNR of an event is relatively lower when observed by ET alone. The rates of acceptance drop to 35.7%(C) / 43.5%(D).
In addition, the statistical characteristics (µ, σ) of SNR for the four catalogs are: A (22.45, 13. With these ready-to-use data, we then evaluate the accuracy of distance measurement.

Error analysis
We uniformly divide the redshift range z ∈ (0, 5) for A and B or z ∈ (0, 3.5) for C and D into 10 bins, and calculate the average fractional error of d L in each bin. The results are plotted as functions of redshift in Figure 3. Again, it is obvious that smaller σ v and more detectors both lead to smaller errors at given redshifts. Specifically, the improvement of accuracy due to more detectors is evident in the whole redshift range, almost halving the uncertainty in each bin. While the impact of smaller σ v turns out to be more evident at low redshifts. This tendency can be explained by the simplified expression eq. (25). At low redshifts, the average SNR is relatively high; thus, the error induced by σ v overwhelms the 1/ρ term, making the alteration of σ v more significant.
Besides, by means of simulation, [35] found that without the help of the EM counterpart, the average relative error of d L in z ∈ (0, 1) was 0.48 for ET and 0.33 for the ET + LISA + B-DECIGO network, whereas, in the presence of SGRB, we have σ ∆d L /d L < 0.1 in the same range for each configuration. This is direct evidence of the benefit of the EM counterpart.
With dashed curve, we also plot σ ∆d L ,0 /d L in the same figure for comparison. In each case, adopting σ ∆d L ,0 causes overestimation, meaning that it is a rather conservative choice for the error analysis.
Another factor usually accounted for is the systematic error due to weak lensing [58,59], which is a major source of error on d L for high-redshift standard sirens Thus, the total error can be expressed as σ ∆d L ,tot = σ 2 ∆d L ,lens + σ 2 ∆d L .   We also plot the lensing errors in Figure 3. In each panel, the contribution of lensing is comparable to the instrumental noises. So far, we have obtained four groups of {z, d L , σ ∆d L ,tot } values. However, these (z, d L ) pairs strictly satisfy the functional relationship eq. (26), which violates the stochastic nature of the universe. Thus, we will generate the luminosity distances by Gaussian distribution N(d L , σ ∆d L ,tot ), which is a good approximation when the error is not too large. The data points as well as associated errors for the four configurations are shown in Figure 4.

ΛCDM cosmology
To determine how many events are required to achieve acceptable precision for cosmological parameters, again we make use of the FIM. Under flat ΛCDM model, the FIM of parameters p = (H 0 , Ω m ) is as follows: where the partial derivatives are calculated according to eq. (26), and the covariant matrix is Σ ΛCDM = Γ ΛCDM −1 .
The errors of parameters are presented in Figure 5. Both σ H 0 and σ Ω m scale roughly as 1/ √ N. For each parameter, the improvement in accuracy from large σ v to small σ v is significant , while for H 0 , the influence of detector number seems moderate at given N. As a matter of fact, the detecting capability of ET is poorer than that of the ET + CE network. Judging from the rates of acceptance, the number of events recognized by ET is only about half of ET + CE in the same detecting period. Thus, the values of N should be different when it comes to the comparison between single and multiple detectors (e.g., N = 2 for ET and N = 4 for ET + CE, with σ v fixed), and it turns out that the impact of detector number is even more significant than that of σ v value. Analyses in the preceding sections are not affected by this factor since only the individual and average errors are involved there.
To achieve 1% accuracy of H 0 , the numbers of events required for the four configurations are about 500, 200, 600, and 300, respectively. With all 1,000 events, H 0 can be constrained to 0.66%, 0.37%, 0.76%, and 0.49%, while the errors of Ω m are 0.010, 0.006, 0.013, and 0.010, respectively. These results are comparable to those of [19], which declared that a subpercent level accuracy on H 0 could be reached by ET.
As supplement, we also present the results of full Bayesian analysis, performed with the public code CosmoMC [60]. Figure 6 and Table 2 show the posterior distributions of parameters p = (H 0 , Ω m ) constrained from four catalogs of 1,000 events. The blue contours correspond to the results obtained by adopting σ ∆d L ,0 . After simple calculation, it is easy to verify that the results of Markov chain Monte Carlo (MCMC) is compatible to those of FIM. And very reasonably, using σ ∆d L ,0 leads to overestimation, regardless of the number of detectors or the value of σ v , consistent with our analysis in sect. 3.2.

Cosmology with generic H(z)
Despite model-independent approaches, GP [61], as a nonparametrization method, can extract information from observational data without any hypothetic functional form. Therefore, GP provides a means to study how accurately a generic function H(z) can be evaluated from data. To this end, we first calculate d L (z) as a GP where µ(z) and k(z, z ) are the mean value and covariance function, respectively. Then, H(z) can be obtained via where " " denotes the derivative of redshift. The only assumption behind eq. (34) is that the geometry of the universe is depicted by the flat Friedmann-Robertson-Walker metric.
Using the open-source code Gaussian Processes in Python (GAPP), we reconstruct H(z) from the four catalogs, and the resulting mean values and errors up to 2σ are presented in Figure 7. Future GWSS has the advantage of deeper detection depth than conventional standard candles, such as SN Ia. For example, the recently released Pantheon sample [62] consists of 1,048 SNe Ia with maximum redshift 2.26, and most of them are within the range of z < 1. While, as is shown in Figure 7, GWSS can constrain H(z) to considerable accuracy at relatively high redshifts. The 1σ precision of 10% can be maintained to redshift 3.12, 3.34, 2.37, and 2.53 by configurations A, B, C, and D, respectively. Besides, the H(z) functions reconstructed from A and B agree with the fiducial model (blue curves), even when z → 3.5, whereas, for C and D, ΛCDM exceeds the 1σ ranges at z = 1.90 and z = 2.15. In conclusion, adopting multiple detectors and tighter prior both contribute to a precise and unbiased measurement. The constraining ability of future GW data can be better demonstrated when compared with current observations. Compiled by [63] are the Hubble parameters measured at different redshifts within z ∈ (0, 2) (observational Hubble data, OHD) via methods such as the cosmic chronometer and baryon acoustic oscillation. The OHD with error bars and the results of GP are shown in the upper panel of Figure 8. For comparison, the mean values and fractional uncertainties of H(z) from GW data and OHD are summarized in the upper and lower panels, respectively. In almost the whole redshift range, H(z) can be measured more accurately from each GW catalog than from OHD.

Further comparisons
The entire pipeline of simulation and data processing is based on several models and hypotheses, such as the distribution of inclination angle and the SFR model. In this section, we make alterations to them and investigate what difference will be brought about.

Variedv vs. fixedv
In this section, we investigate the impact of varied mean value of v in the prior, i.e., replacing eq. (23) with v being the "true" source parameter. This analysis is performed to simulate a situation wherev is already determined from the SGRB observation with precision σ v ; thus, the data of SGRB and GW are employed in a more joint manner. In the simplified regime where there is only one stationary interferometer, the consequence of using eq. (35) can be derived analytically: By numerical calculation, it is easy to confirm that K > 1 for the typical values of v 1 , v 2 , and ρ; thus, varyingv usually increases error compared with eq. (25). It should be noted that eq. (25) is estimated in the face-on limit, while eq. (36) is valid whether v is 1 or not; thus, the performances of these two priors remain to be compared through simulation. Besides, in practice, the uncertainties of v constrained from SGRB data vary among events; thus, using a single value σ v is a rather rough choice. Following the prescription of error budget given in sect. 3.2, we apply eq. (35) to the simulated GW catalogs, and the instrumental errors are shown in Figure 9 as solid curves. For comparison, we also plot the results of eq. (23) in the same figure with dotted curves. It is obvious that varyinḡ v increases the errors of d L . Besides, for the cosmological parameters p = {H 0 , Ω m }, we have σ H 0 /H 0 = 0.74%, 0.37%, 0.84%, 0.51% and σ Ω m = 0.011, 0.006, 0.014, 0.010, respectively, which are also slightly larger than those of fixedv. Moreover, the influence of adopting eq. (36) is more evident when σ v = 0.05 since v is allowed to vary in a wider range in this case.

Other star formation rate models
The simulation of GW sources is based on the hypothesis that the NS formation follows the SFR. While, as a matter of fact, the high-redshift SFR is uncertain at present. Pioneering studies [64][65][66] have revealed that GRBs observed by Swift provide a biased measurement of the SFR history, with an enhancement of ∼ (1 + z) 0.5 . [64] reported that this factor could be readily explained if GRBs occur primarily in lowmetallicity galaxies proportionally more numerous at earlier times. However, after considering this trend, the SFR beyond z = 4 derived from GRBs is still much higher than that from other surveys, such as ultraviolet (UV) and far-infrared (FIR) measurements. [65] provided an explanation for the highredshift GRB rate excess by considering the GRBs produced by rapidly rotating metal-poor stars from low masses.
It is beyond the scope of this paper to explore the physical origins of the discrepancies between different SFR surveys. On the contrary, apart from Fiducial + PopIII (eq. (29)), we adopt two additional SFR models, one inferred from GRBs [66] (dubbed Kistler): The simulated GW catalogs are compared in Figure 10. For clarity, we only present the events obtained under configuration B. The red, blue, and black histograms represent the distributions of GW sources generated from Hopkins, Kistler, and Fiducial + PopIII, respectively, and curves with the same colors correspond to the theoretical BNS merger rates. For both Hopkins and Kistler, the peaks of P(z) come at lower redshifts (z = 1.1) than Fiducial + PopIII (z = 1.5), and the high-redshift merger rates are slightly larger, especially for Kistler, due to the so-called "GRB rate excess" phenomenon.
Under configuration B, the uncertainties of cosmological parameters are σ H 0 /H 0 = 0.295%, 0.349%, 0.367% and σ Ω m = 0.0057, 0.0064, 0.0064, for Hopkins, Kistler, and Fiducial + PopIII, respectively. The errors of each parameter are in the same order of magnitude. The comparison between simulations agrees with intuition, i.e., parameters can be better constrained from catalogs with more low-redshift (high-SNR) events.

Hopkins Kistler
Fiducial + PopIII Figure 10 The redshift distributions of gravitational wave events (configuration B) simulated according to different star formation rate models. The red, blue, and black histograms represent the results of SFR models, Hopkins, Kistler, and Fiducial + PopIII, respectively, and the curves with the same colors correspond to the theoretical BNS merger rates.

Other short gamma-ray burst models
In former sections, we have assumed that SGRB follows a Gaussian distribution peaked at v = 1, in concordance with [21,40]. The modeling of SGRB plays a crucial role in the process of parameter estimation, and several models, such as the jet and cocoon scenarios [46], have been proposed to describe the SGRB features in the afterglow of GW. Therefore, we discuss this issue by comparing our model with the Gaussian structured jet profile which is a Gaussian distribution on ι with standard deviation σ ι . [46] constrained σ ι = 0.091 +0.037 −0.040 rad, 0.057 +0.025 −0.023 rad, and 0.076 +0.026 −0.027 rad from different data combinations, whereas [47] reported σ ι = 4.7 +1.1 −1.1 deg. These values are very close to our choice of σ v = 0.003; thus, we adopt σ ι = cos −1 (1 − 0.003) = 0.077, replace eq. (23) by eq. (39), and repeat simulations B and D.  Figure 11 Relative errors of d L as functions of redshift. Results obtained by adopting p ι (ι) (or p v (v)) are plotted with solid (or dotted) curves, and blue (or red) curves correspond to configuration B (or D).
It appears that the instrumental errors σ ∆d L /d L (z) of different priors are indistinguishable ( Figure 11). Thus, we further calculate the uncertainties of ΛCDM parameters. With eq. (39), we obtain σ H 0 /H 0 = 0.34%, 0.58% and σ Ω m = 0.006, 0.011 for B, D, respectively. The results of adopting p ι (ι) is only barely different from those of p v (v) since they both indicate that ι is confined in a small range (< 5 deg) around 0.

Conclusion
In this paper, we address the problem of degeneracy between luminosity distance and inclination angle in the estimation of GW parameters, especially near the face-on limit. For GW-SGRB events, the sky locations of sources can be pinpointed with techniques such as identifying the host galaxies, and the mass parameters can be determined with exquisite accuracy from the phase of GW; thus, d L and v = cos ι are treated as independent of other parameters. We first calculate the FIM of these two parameters and then consider a Gaussian prior on the inclination based on assumed knowledge about SGRB. This pipeline is tested via simulation. Four catalogs, each with 1,000 events detected by 3G GW detectors, are simulated, and the settings are {A: ET + CE / σ v = 0.05, B: ET + CE / σ v = 0.003, C: ET / σ v = 0.05, D: ET / σ v = 0.003}. Although future SGRB observations may give σ v other than these two values, our research may give a clue to the changes in results when σ v varies.
It turns out that the detector network, ET + CE, can recognize more and farther events than a single detector ET in the same detecting period. By analyzing the error of d L , we find that using more detectors and tighter prior on the inclination can both improve the precison of measurement. Also considered is the performance of a widely adopted formula σ ∆d L ,0 = 2d L /ρ, which overestimates the error of d L for each configuration.
The simulated catalogs are further applied to constrain cosmological parameters. For ΛCDM cosmology, 500, 200, 600, and 300 events are required for the four configurations to achieve 1% H 0 accuracy. With all 1,000 events, H 0 is constrained to 0.66%, 0.37%, 0.76%, and 0.49%, whereas the errors of Ω m are 0.010, 0.006, 0.013, and 0.010, respectively. The MCMC method produces results compatible with those of FIM. Regarding the Hubble parameter as a free function of redshift and using the GP method, we find that H(z) can be measured with uncertainties less than 10% up to redshifts 3.12, 3.34, 2.37, and 2.53, indicating that future GWSS can probe cosmology at higher redshift than the observed standard candles to date. The advantage of future GWSS over current OHD is also obvious. Besides, adopting more detectors and tighter prior both contribute to a precise and unbiased result.
In addition, we have also considered the effect of altering some of the formulae and models used in the simulation and data processing, such as SFR and the distribution of inclination angle. Detailed analyses can be found in sect. 4.
In this paper, the prior distributions of inclination are based only on limited observations. With upcoming detections on SGRB and GW, as well as more reliable estimations of the joint event rate, we expect that more rigorous results can be achieved.