On Dark Gravitational Wave Standard Sirens as Cosmological Inference and Forecas

Ju Chen,Changshuo Yan,Youjun Lu,Yuetong Zhao,and Junqiang Ge

1 National Astronomical Observatories,Chinese Academy of Sciences,20A Datun Road,Beijing 100101,China; yancs@nao.cas.cn,luyj@nao.cas.cn

2 University of Chinese Academy of Sciences,19A Yuquan Road,Beijing 100049,China

Abstract Gravitational wave (GW) signals from compact binary coalescences can be used as standard sirens to constrain cosmological parameters if their redshift can be measured independently by electromagnetic signals.However,mergers of stellar binary black holes (BBHs) may not have electromagnetic counterparts and thus have no direct redshift measurements.These dark sirens may be still used to statistically constrain cosmological parameters by combining their GW measured luminosity distances and localization with deep redshift surveys of galaxies around it.We investigate this dark siren method to constrain cosmological parameters in detail by using mock BBH and galaxy samples.We find that the Hubble constant can be constrained well with an accuracy ≾1%with a few tens or more of BBH mergers at redshift up to 1 if GW observations can provide accurate estimates of their luminosity distance (with relative error of ≾0.01) and localization (≾0.1 deg2),though the constraint may be significantly biased if the luminosity distance and localization errors are larger.We also introduce a simple method to correct this bias and find it is valid when the luminosity distance and localization errors are modestly large.We further generate mock BBH samples,according to current constraints on BBH merger rate and the distributions of BBH properties,and find that the Deci-hertz Observatory (DO) in a half year observation period may detect about one hundred BBHs with signal-to-noise ratio ϱ ≿30,relative luminosity distance error ≾0.02 and localization error≾0.01 deg2.By applying the dark standard siren method,we find that the Hubble constant can be constrained to the～0.1%–1% level using these DO BBHs,an accuracy comparable to the constraints obtained by using electromagnetic observations in the near future,thus it may provide insight into the Hubble tension.We also demonstrate that the constraint on the Hubble constant applying this dark siren method is robust and does not depend on the choice of the prior for the properties of BBH host galaxies.

Key words: gravitational waves–cosmological parameters–cosmology

1.Introduction

The Hubble constant (H0),which measures the current expansion rate of the Universe,plays a fundamental role in state of the art cosmological models.An accurate measurement of H0is crucial for understanding the Universe,and numerous methods have been developed to achieve this goal(Baade 1944;Blandford&Narayan 1992;Riess et al.1998;Perlmutter et al.1999; Freedman et al.2001; Spergel et al.2007; Freedman &Madore 2010; Beutler et al.2011; Jackson 2015).Most notably,with the latest Planck cosmic microwave background(CMB) data,H0is determined to be 67.4±0.5 km s-1Mpc-1with a precision of～1%,assuming a standard ΛCDM cosmology (Planck Collaboration et al.2018),while it is measured to be 74.03±1.42 km s-1Mpc-1by relying on the distance ladder measurement with a precision of～2% (Riess et al.2018,2019).These estimates are discrepant by more than 4σ.The result measured with lensed quasars recently is also different from the Planck measurement by 3σ,and together with distance ladder,it suggests a 5σ tension between H0measurements obtained by early-and late-Universe probes(Wong et al.2020).This tension may be an indicator of unknown systematic bias in current methods or new physics beyond our current standard cosmological models.It is therefore extremely important to find new and independent ways to measure the Hubble constant.

Compact binary coalesces may be used as “standard sirens” to measure H0and other cosmological parameters (Schutz 1986),if their redshift(information)can be obtained by independent method(s),as the luminosity distances of these sources can be directly measured from their gravitational wave(GW)signals.This method has been intensively investigated in the literature by using mock data(Holz&Hughes 2005;Nissanke et al.2010;Zhao et al.2011;Taylor & Gair 2012; Nissanke et al.2013; Tamanini et al.2016;Liao et al.2017;Chen et al.2018;Li et al.2019;Zhao et al.2020;You et al.2021;Wang et al.2021).The advantage of this method has been demonstrated by using observations of the first binary neutron star(BNS)merger GW170817,the redshift of which was measured via the discovery of its electromagnetic (EM) counterpart and host galaxy(NGC 4993)(Abbott et al.2017;Hotokezaka et al.2019).

To apply the“standard siren”method,it is necessary to have redshift measurements of binary coalescences by utilizing their EM counterparts since they cannot be obtained directly from the GW observations due to the degeneracy between chirp mass and redshift.However,it is not an easy task to identify their EM counterparts and/or host galaxies due to large uncertainties in the sky localization given by the GW signals and the small opening angle/faintness of the EM emission.Moreover,a large fraction of the GW sources detected by ground-based detectors would be mergers of binary black holes(BBHs),which are even not expected to be accompanied by bright EM counterparts.

The planned third generation ground-based GW detectors(e.g.,Einstein Telescope: ET,and Cosmic Explorer: CE) and future space GW detectors (e.g.,DECIGO and Deci-hertz Observatory,DO) are expected to detect a large number of compact binary mergers with high signal-to-noise ratio (S/N)at redshift up to 10,and will yield～104～108detections of BBH events per year (Abernathy et al.2011; Abbott et al.2017;Arca Sedda et al.2020).Most of these mergers,if not all,will not have detectable EM counterparts and thus redshift information.It is then of great interests to investigate whether these dark “standard sirens” can be used to probe cosmology and constrain cosmological parameters.

Schutz (1986) originally considered that each galaxy(cluster) within the localization error volume can be taken as a potential host candidate,and by combining different GW events,only one galaxy(cluster)in each error volume will give consistent constraints on cosmological parameters.The effectiveness of this method has been discussed recently in various cases including both ground-based (Del Pozzo 2012; Chen et al.2018; Gray et al.2020) and space-based (MacLeod &Hogan 2008;Petiteau et al.2011;Del Pozzo et al.2018;Wang et al.2020) detectors.Yu et al.(2020) further investigated the capabilities of identifying host galaxy groups instead of host galaxies in considering the hardness of carrying out deep galaxy survey.Such a method is also applied to the real observational data as presented in Fishbach et al.(2019) for GW170817 without considering the EM counterpart,and The DES Collaboration et al.(2019) for GW170814.The LIGO Scientific Collaboration et al.(2021) also presented the result obtained by combining observations of the current standard siren GW170817 and all other detected dark GW events.

In this paper,we investigate the validity and limitation of relying on dark “standard sirens” (mainly mergers of stellar BBHs in this paper) to constrain cosmological parameters and quantify the accuracy of such a method and the selection criteria for the sample of dark “standard sirens” to be used for cosmological application.The paper is organized as follows.In Section 2,we introduce the framework used in our analysis.In Section 3,we describe the detailed method for generating GW sources and host galaxy mock data based on a simple model for the evolution merger rate density and the Millennium database.In Section 4,we check the robustness of the method,and discuss the systemic bias we noticed in detail and also introduce a correction factor.In Section 5,we further investigate whether the DO could provide a large number of GW events with sufficiently smaller errors of luminosity distance measurements and localization and test the ability of DO observations to constrain the cosmological parameters via dark sirens.Discussions are given in Sections 4.4 and 6.Conclusions are summarized in Section 6.

Throughout the paper,we assume a flat ΛCDM universe,with the Hubble constant H0,matter fraction of Ωmand cosmological constant (1-Ωm).We adopt the conventional dimensionless parameter h=H0/100 km s-1Mpc-1to denote the Hubble constant.

2.The General Framework

In this section,we briefly introduce the general framework to infer cosmological parameters by utilizing the dark standard sirens,e.g.,mergers of BBHs without EM counterparts and even mergers of faraway BNSs with too faint EM counterparts to be detected.GW signals provide estimates of the luminosity distances and sky localization of these dark sources but not the redshift.EM galaxy surveys can provide host candidates for these GW sources in the cosmic volume defined by the errors of luminosity distance and sky location of each source.The combination of the GW measurements and galaxy surveys may enable constraints on cosmological parameters by using these dark standard sirens (Del Pozzo et al.2018,see also Del Pozzo 2012; Petiteau et al.2011).

For simplicity,we assume the sky localization area for a GW event is circular in the following analysis although the one given by GW observation is normally elongated and even irregular.Considering that our results are only affected by the distribution of galaxies along the redshift direction but not the direction transverse to the line of sight (LOS),the approximation of a circular localization area should not have any significant effect on our results,though one must consider this shape in practical studies.

With the information from both GW observations and galaxy surveys described above,the posterior probability distributions for the cosmological parameters,including the Hubble constant(H0or h),fraction of matter content (Ωm) and others,may be estimated by

according to the Bayes’ theorem.Here evidence D=(D1,…,Di,…,Dn) with Direpresenting the observational data,including the data set digiven by GW observations and the redshift zi=(zi,1,…,zi,j,…)of host galaxies in the localization error volume of each GW source i,I signifies all the additional information and assumptions,p(H|I) is the prior probability distribution for cosmological parameters,p(D|H,I) is the likelihood function for evidence D,and p(D|I)is the probability distribution for evidence D under all possible parameter spaces and can be treated as a normalization constant here.Since GW events should be independent of each other,the likelihood above can be rewritten as the product of the likelihoods for individual events Di,i.e.,

The likelihood for a single GW event is obtained by marginalization over all different parameters x of the source system,i.e.,

For the consideration in the present paper,these parameters include the celestial coordinates θ,φ,luminosity distance dLand redshift z,i.e.,x ≡{θ,φ,dL,z}.We assume that the distribution of GW source luminosity distance dLis independent of sky location (θ,φ),then we have

The likelihood for each GW event p(Di|x,H,I)is assumed to be Gaussian distribution centered on the luminosity distancewith a standard deviation ofσdLiinferred from the GW signal,i.e.,

Here the luminosity distance of the host galaxies in the localization error volume of each GW source i is expressed as

With the assumption that GW sources are hosted in galaxies,the joint distribution of other parameters p(θ,φ,z|H,I) can be obtained by summing up the probability of all galaxies located in the localization error volume being the host.Here we simply choose a sum of δ-functions over all individual galaxies,i.e.,

where j goes through all the cataloged galaxies located in the localization error volume of each GW event Di,m is the total number of these galaxies and it may be significantly different for different GW events,pjis the weight assigned to each of these galaxies according to additional information,if available,about the probability distribution of host galaxy properties,and we assume it does not depend on θj,φj,or zj.In principle,pjis not known as a prior.The simplest way to consider this probability is to assume that the BBH merger rate for each galaxy is proportional to its stellar mass M*.In this case,the probability of a galaxy with M*to be the host is

where dN/dM*is the stellar mass function of galaxies,which in our case for all galaxies in the error volume of each GW event is ∑jδ(M*-M*,j).

Note here we omit the additional subscript i for each GW event in the above equation and hereafter for simplicity,but we should keep in mind that for each GW event Di,m can be quite different.

A more realistic distribution of host galaxy properties may be statistically inferred by future GW detection of a large number of BBH mergers.Such a distribution can also be derived by considering detailed formation and evolution of BBHs along with galaxies in the Universe (Cao et al.2018;Artale et al.2019).In this paper,we simply use pj∝M*,j to assign the host galaxy for each GW event when generating the mock sample (for details see Section 3).For the inference of cosmological parameters,we also assume pjis known as a prior,and set the weight according to galaxy stellar mass if not otherwise stated.The effects of different choices of pjon the cosmological parameter inference will be further discussed in Section 4.4.

According to Equations (4) to (8),the likelihood for an individual GW event is then

By integrating over parameters dL,θjand φj,we have

In the above Equations (8) and (11),the uncertainty in the redshift measurement for each galaxy zj,e.g.,due to peculiar velocity,is not considered.However the influence of peculiar velocities is not negligible for nearby galaxies.To account for the effect of peculiar velocity,we assume the uncertainty of redshift introduced by peculiar velocity is a normal distribution around the true redshift,then we have

wherezj,obsis the observed redshift including peculiar velocity and can be obtained from the mock catalog (Henriques et al.2012).The standard deviationσzpvis set to 0.0015 for low redshift galaxies (z＜0.1) and 0.0015+0.0005z for galaxies with higher redshift.These values are chosen based on the distribution of redshift measurement error for all galaxies in the mocked galaxy catalog in Henriques et al.(2012).

With the above assumption about redshift distribution of galaxies,the likelihood for a GW event becomes

The upper and lower limits (zmin,zmax) of the integration are estimated as described in Section 3.3 for the selection of host galaxy candidates.

To investigate the prospects of the dark siren method described above for future GW observations,below we generate mock samples of GW events with luminosity distance and localization measurements.Correspondingly,for each GW event we can get a list of galaxies as host candidates from the mocked galaxy catalog given by Henriques et al.(2012).Then,we adopt the above framework for Bayes’s inference to reconstruct cosmological parameters from the mock data.We adopt the nested sampling algorithm provided by dynesty(Speagle 2020)to carry out the posterior distribution sampling.We first investigate the dark siren method by using mock data with a set of assumed uncertainties for the luminosity distance measurements and localization via GW observations,and then test with mock data simulated for future DO observations(Arca Sedda et al.2020).

To generate mock catalogs of GW events,we assume a flat ΛCDM Universe with h=0.73,Ωm=0.25.The parameters are chosen to be the same as the Millennium Simulation for consistency as the host galaxies are selected by using the galaxy catalogs obtained from the Millennium Simulation.

3.Mock Data for GW Events and their Host Galaxies

3.1.Stellar Binary Black Holes

We first generate mock samples of stellar BBHs at z＜3 by adopting the BBH merger rate density (and its evolution)inferred from the LIGO/Virgo observations after O3a (The LIGO Scientific Collaboration et al.2020),i.e.,

where R0is the local merger rate density.According to the power law plus a peak model in The LIGO Scientific Collaboration et al.(2020),we have R0=19.1 Gpc-3yr-1and k=1.3.We ignore the difference of the cosmological model adopted for the estimation of R0from the “true”cosmological model (Millennium Simulation) considered in this paper for mocking GW events and their host galaxies,as current estimation on R0still has a large uncertainty.We assume that the BBH merger rate density follows Equation(14)up to z=3 although it is obtained by only using LIGO/Virgo detected BBHs at low redshift.

3.2.Host Galaxies of GW Events

In this subsection,we describe how to generate mock samples of GW events and their host galaxies.According to Equation (14),we can randomly generate a given number of BBH mergers with assigned redshiftz′ .We adopt the virtual observed light cones of Henriques et al.(2012),constructed from the semi-analytic galaxy formation model (Guo et al.2011)based on the Millennium Simulation(Springel 2005).In these light cones,galaxies with different stellar masses (or other properties) may not have the same probability to be the host galaxy of a GW event.In order to account for the dependence of a GW event on galaxy stellar mass,as discussed in Section 2,we generated a catalog of GW host galaxies from the simulated one given by Henriques et al.(2012),by convolution of the probability of a galaxy with mass M*to be the GW host p(M*).By default,we will choose the distribution simply as p(M*)∝M*,and we will further discuss the influence of p(M*) in Section 4.4.For each BBH merger atz′,we may then find a galaxy in the new catalog for GW hosts with redshift z that is the closest one to the BBH merger atz′ .We replace the redshift of the GW eventz′ by z,and its host has a mass of M*.By doing so,we can generate catalogs of GW events with known host galaxies and perform further analysis.According to the light cones (assuming periodical boundary condition),we can get the sky location,stellar mass,real redshift z,observed redshift zobs(including the effect of peculiar velocity),and apparent magnitude of each galaxy for each GW event and its host galaxy.Note here that galaxies with M*＜109M⊙are excluded from our consideration,since they are rare and may be biased due to the resolution limitation of the Millennium Run.

As a test of the method,we first generate some mock samples(both at redshift z＜1 and z＜0.1)of GW events with assumed fixed uncertainties of the luminosity distance measurements and localization in Section 4.After that,we investigate the prospective of future observation with mocked detection of DO in Section 5.In these cases,we first obtain the luminosity distanced′Lfrom the true redshift of the host galaxy of each GW event under the cosmological model (h=0.73,Ωm=0.25).We then get the observed luminosity distancefor these GW events by adding a random error tod′L,and the errors are chosen to be Gaussian distribution with standard deviationσdLgiven by mock observations.

3.3.Selection of Host-galaxy Candidates for GW Events

Given the luminosity distance and its uncertainty±δdLof a GW event,we can get the corresponding redshift range for any given set of cosmology parameters.If the cosmology parameters are not known prior,the error range for the redshift of the GW event would be significantly larger as one must account for the uncertainties in the cosmological parameters.

We choose the prior for h and Ωmas flat distributions in[0.6,0.8] and [0.01,0.5],respectively.With a given set of (h,Ωm),according to the flat ΛCDM model,the luminosity distance and its error±σdLof GW event can be converted into a redshift range[z-,z+].In this paper,we adopt a3σ dLrange to determine the redshift error range.Under the assumed priors for h and Ω,the minimum value of z-is from the model with (h,Ωm)=(0.6,0.01),and the maximum value of z+is from the model with (h,Ωm)=(0.8,0.5).Furthermore,galaxy surveys normally give the observed redshift of each galaxy,which deviates from the true one because of the contamination by peculiar velocity.To account for this contamination,the redshift error range is enlarged from (z-,z+) to(zmin,zmax),which includes the errors induced by the peculiar velocity aszmin=z--3σ zpv,zmax=z-+3σ zpv,withσzpvalso set to 0.0015 for low redshift galaxies (z＜0.1) and 0.0015+0.0005z for galaxies with higher redshift.

With the redshift range determined above(without fixing the cosmological parameters h and Ωm),the error volume can be obtained by considering both the localization error ΔΩ and redshift rangezmin,zmax[ ].All galaxies falling within the error volume are then selected as potential hosts.

For the galaxy catalogs given by Henriques et al.(2012),we use the all-sky galaxy map for GW events at low redshift(z＜0.15),while the pencil-beam light cones assume a periodic expansion for events at high redshift.We also set a threshold on the i-band apparent magnitude limit of host galaxies as 25 at high redshift.

4.Robustness Check of the Method

4.1.Systematic Bias in the Method

We first investigate the dependence of cosmological constraints on the measurement errors of luminosity distance and localization of GW events.Note,for real observation the errors are normally smaller for nearby events and larger for more distant events.Here to demonstrate the dependence of the method on the observation errors,we will first ignore this and simply assume that all GW events have the same uncertainties in luminosity distances and localization.More realistic mock observation will be presented in Section 5 for DO.

Figure 1 displays the posterior probability distributions of cosmological parameters h and Ωmobtained by using 50 mock GW events at redshift z ≤0.1 (left panel) and z ≤1 (right panel).The uncertainties(90%level)of the luminosity distance measurements from the GW signals for all events are set as δdL/dL=0.01,while the sky localization uncertainties (90%level) are set as ΔΩ=5 deg2and 0.1 deg2for GW events applied in the left and right panels,respectively.As shown with the blue color results,h can be constrained with a precision of～1% even if only 50 dark sirens with redshift up to 1 are adopted.However,Ωmcannot be constrained well,especially if only low redshift (z＜0.1) GW events are available,which is consistent with Del Pozzo et al.(2018).Adopting 50 GW events with redshift up to 1,the constraint on Ωmcan be improved but only reach a precision of～8%.

Figure 1.Posterior probability distributions of h and Ωm obtained from mock GW events at redshift z ≤0.1 (left panel) and z ≤1 (right panel).The relative uncertainty for the luminosity distance δdL/dL is set to 0.01,and the sky localization uncertainty ΔΩ is fixed at 5 deg2 for mock GW events adopted at z ≤0.1 and 0.1 deg2 at z ≤1.The blue and red contours in the bottom-left sub-figure of each panel signify the 1,2 and 3σ level constraints obtained from the method described in Section 2 without and with consideration of the correction factor introduced in following,respectively.The blue and red curves in the top and bottom-left sub-figures of each panel correspondingly depict the marginalized one-dimensional probability distributions of h and Ωm(in the horizontal direction),respectively.Vertical and horizontal dashed lines indicate the input values of h and Ωm,and the median values of recovered parameters are shown as“+”symbols.This figure affirms that the cosmological parameters can be constrained well when the luminosity distances and localization of the GW sources are measured with high precision,and with the correction factor α(H) introduced in Section 4.2,the constraints are only slightly improved.

Figure 2.Legend is similar to Figure 1,except that the measurement errors of the luminosity distances and localization areas for GW sources are set to (δdL/dL,ΔΩ)=(0.1,5 deg2)and(0.01,1 deg2),respectively.This figure demonstrates that the inferred cosmological parameters are highly biased when the measurements of the luminosity distances and localization of GW sources are not precise enough (see Figure 1),and the constraints can be improved significantly once the bias correction discussed in Section 4.2 is considered.

Figure 3.Recovered cosmological parameters (median value case of 5 runs) from 50 mock GW events at z＜0.1 by assuming different luminosity distance and localization measurement errors.In the left panel,the localization errors are all assumed to be 1 deg2,but the relative luminosity distance errors are assumed to be δdL/d,L=0.01,0.05,0.1,0.15 and 0.2 as indicated,respectively.In the right panel,the relative luminosity distance errors are all assumed to be 0.01,while the localization errors are assumed to be 0.01 deg2,0.1 deg2,0.1 deg2,1 deg2,5 deg2 and 10 deg2,as indicated,respectively.The error bars associated with each symbol mark the 1σ error range of the constraint.Blue and red stars indicate the results obtained without and with consideration of the bias correction(see Section 4.2),respectively.This figure suggests that the constraints can be improved significantly once the bias correction is considered for those samples with relatively modestly large luminosity distance and localization measurement errors.

Figure 2 displays the constraints obtained from 50 mock GW events with(δdL/dL,ΔΩ)=(0.1,5 deg2)at z ≤0.1(left panel)and (0.01,1 deg2) at z ≤1 (right panel).By comparing Figures 2 and 1,we can see that the constraint on h becomes not only much less tight but also biased away from the input value,if the measurement errors of dLand/or ΔΩ are large.

Figure 3 shows the dependence of the precision of h and Ωmconstraints on the luminosity distance and localization measurement errors of the GW events(totally 50).As signified with blue markers,it is clear that the larger the luminosity distance errors (left panel) or the larger the localization errors(right panel),the poorer the constraint on h,and the constrained value of h is even significantly biased away from the input value when δdL/dL≿0.1 and ΔΩ ≿1 deg2.

To investigate the source for the systematic bias revealed in Figures 2 and 3,we check the likelihood for a single GW event

We arbitrarily set some galaxies uniformly distributed in the redshift error range [z-,z+] of a GW event and assume that each one has an equal probability to be the GW host.Then we calculate the likelihood according to the above equation for different cosmology parameters.Figure 4 displays the results for such a likelihood for cosmological models with various h but fixed Ωm(0.25).In this figure,each gray line indicates the probability distribution of parameter h given by a possible GW host galaxy with an assigned redshift in the redshift error range zj,and the solid line represents the summation of the likelihood of all individual possible host galaxies (signifying the total likelihood for single event) with the peak rescaled to 1 (see Equation(15)).Apparently,the total likelihood peaks at a value of h larger than the input one and the likelihood distribution is not symmetric.The main reason is that the dependence of dLon cosmological parameters (h,Ωm) is nonlinear and the likelihood for individual host candidates is wider at large h than that at small h,which causes the obtained constraints on the values of h and Ωmto deviate from the input ones.

Figure 4.Likelihood of individual GW events with different h but fixed Ωm with uniform distribution of galaxies in the redshift error range.The dotted lines show the contributions from each galaxy in the redshift error range by fixing Ωm=0.25,and the solid line marks the summation of the contributions from all these galaxies with peak rescaled to 1.The vertical dashed line indicates the input value of the Hubble parameter.The dash–dotted line is the likelihood with consideration of the bias correction (see Section 4.2).

Figure 5.Galaxy distribution in the error volume/redshift error range of an example GW event (left panels) and the corresponding likelihood (right panels;Ωm=0.25).Top and bottom panels show the two cases with settings of (δdL/dL,ΔΩ)=(0.01,1 deg2) and (0.1,5 deg2),respectively.The vertical dotted line indicates the true host redshift.The vertical solid line signifies the redshift converted from the best-fit value of luminosity distance dL given by the GW observations and the two dashed lines indicate the redshift ranges corresponding to the 3σ range of dL,assuming the input cosmological model.The dash–dotted vertical line indicates the redshift range after accounting for the uncertainty of cosmology parameters and redshift uncertainty described in Section 3.The gray lines show the likelihood contributed by each galaxy (weighted by mass) in the error volume,and the black line is the summation of all galaxies’ contributions.The lines for each galaxy are scaled by their maximum value.The peak value of total likelihood is also scaled to 1.

In reality,galaxies in the redshift error range of a GW event may be clustered rather than uniformly distributed as assumed above.Figure 5 features the distribution of galaxies in the error redshift range of a typical GW event,which is clearly showing some cluster structures.One may also note that the true host may be not necessarily located in one of the clusters and at the center of redshift error range.When the luminosity distance error is large and thus error volume is large (bottom panels),many galaxies are potential hosts.The likelihoods of these galaxies are wide and they overlap each other,which leads to a flat total likelihood over a large h range with the peak biased to large h.When the luminosity distance error is small and thus the error volume is small (top panels),the number of galaxies that can be the potential host is limited.The likelihoods of these galaxies are narrow and the total likelihood of a GW event has several sharp peaks.By combining the total likelihood from different GW events,the contributions from the peaks without a “true” host cancel out and the one that is the “true” host is singled out,and thus a strong constraint on the cosmological parameters,including h and Ωm,can be obtained.

We emphasize here that,in this dark siren method,the total likelihood for a single event is obtained by the summation of the likelihood of h from individual galaxies in the error volume of each GW event.It is affected by the asymmetric nature and change in the width of the likelihood from individual galaxies,which is different from the case of standard sirens,where the corresponding effect is little because the estimation is obtained only by multiplying (rather than summing) the likelihoods contributed from individual host galaxies.A similar bias effect has been reported recently in Laghi et al.(2021)and Muttoni et al.(2021).They found that low S/N events,with larger measurement errors,are not only less informative for inference of the cosmological parameters,but also tend to produce biased results.A non-Gaussian posterior can be helpful to mitigate this bias effect.A more realistic posterior may be considered for real GW events when using them as dark sirens,though this is computationally expensive.In the present work,the likelihood is either assumed to be Gaussian or derived with the Fisher information matrix (FIM)approximation.The FIM approximation is only valid when S/N is high in which case it also assumes asymptotic Gaussian distributions for the errors of those estimated parameters.

4.2.Correction of the Bias

In order to reduce the systematic bias,here we introduce a factor α(H) to re-balance the likelihood

The likelihood for each individual GW event now can be written as

For more realistic cases,we can still assume that galaxies are uniformly distributed in volume as long as the redshift range considered is not large.Therefore,

With the above simple assumption,we re-run the previous cases by including the bias correction factor α(H),and also depict the results in Figures 1,2 and 3.Apparently,the constraints are only slightly improved by including the correction factor for those cases with high precision observations.However,the constraints can be significantly improved for those cases with low precision observations,and the biases can be removed for cases with modest measurement error.While we also notice that with larger observational error (e.g.,(δdL/dL,ΔΩ)=(0.1,1 deg2) for the high redshift case) the results are still biased,especially for Ωm,which need more investigation.

Figure 6.Constraints on h vs.the number of GW events for mock events at z ≤0.1 with (δdL/dL,ΔΩ)=(0.01,1 deg2).Top and bottom panels display the results with and without consideration of the peculiar velocity induced uncertainty in the redshift measurements of galaxies in the mocked catalog respectively.

4.3.Effect of Event Numbers and Peculiar Velocity

One may expect that the constraints from dark sirens,with given measurement errors of the luminosity distance and localization,can be significantly improved if the number of available dark sirens increases to a much larger value.In order to check this,we further investigate the variation of constraints on the cosmological parameters with increasing number of BBHs that are adopted for the inference.

Figure 6 displays the constraint on h obtained from 50,100,200,and 300 mock GW events at z＜0.1 with (δdL/dL,ΔΩ)=(0.01,1 deg2) with (top panel) and without (bottom panel) consideration of the peculiar velocity induced error to galaxy redshift.It is clear that the constraint on h becomes tighter with increasing number of adopted GW events and the statistical error of the constraint shrinks as expected (bottom panel),when both the luminosity distance and localization errors are sufficiently small.However,the best-fit value of h may be significantly biased away from the input value with increasing number of GW events adopted in the analysis when the peculiar velocity induced redshift error is accounted for,though the statistical error indeed becomes smaller(top panel).If ignoring the peculiar velocities of host candidates,i.e.,the redshift error range for a GW event is solely determined by δdL/dL,then the constraint on h is improved with an error shrinking as expected (∝N1/2).This suggests that the constraint by BBHs at z＜0.1 on h may still be able to reach a precision substantially smaller than 1% if the peculiar velocities of host candidates can be corrected on a galaxy-bygalaxy basis using the linear-multi-attractor model similar to that in Freedman et al.(2019).

Figure 7 shows constraints on h obtained from 50,100,200 and 300 mock GW events at z＜1 with(δdL/dL,ΔΩ)=(0.01,0.1 deg2).It is clear that h is constrained well and the precision of the constraint improves with increasing number of GW events adopted.Since the BBHs in these samples span a large redshift range,the effects due to peculiar velocities are negligible.These suggest that only those GW events at high redshift with high precision dLand ΔΩ measurements,e.g.,δdL/dL≾0.01 and ΔΩ ≾0.1 deg2,can be used as dark sirens to tightly constrain the Hubble constant (and other cosmological parameters) with a precision＜1%.

Figure 7.Constraints on h vs.the number of GW events for mock events at z ≤1 with (δdL/dL,ΔΩ)=(0.01,0.1 deg2).

4.4.Effect of the Host Galaxy Distribution Prior

We adopt the above settings for the weight p(M*)to generate mock catalogs of 50 GW events and their host galaxies at redshift z＜0.1.Similar to the previous analysis,we assume δdL/dL=0.01 and ΔΩ=1 deg2.For each mock catalog with a given setting for pj,we use Bayes’ inference to constrain the cosmological parameters by assuming either the “true” p(M*)or the other two types.Our results are listed in Table 1,which suggest that the setting on the prior for the dependence of pj(probability of a galaxy to be the host)on galaxy properties has little effect on the constraints of cosmological parameters (h and Ωm).The reason might be that (1) the redshift error range for each GW event is large and different types of galaxies can always be selected with representative numbers,and (2) the peculiar velocities of galaxies may lead to mixing of the observed redshift distributions for different types of galaxies and thus there will always be different types of host candidates that have redshift close to the “true” host.Therefore,we conclude that the Bayes’inference of cosmological parameters is robust even if no information about the properties of host galaxies is provided.This supports the proposal of using galaxy groups as a surrogate for individual host galaxies to constrain cosmological parameters in Yu et al.(2020).

Table 1 Effects on the Constraints to h and Ωm by Choosing Different Settings for the Properties of GW Host Galaxies

5.Expected Cosmological Constraints from Deci-hertz Observatory

Then,we can generate a mock sample of BBHs by the Monte Carlo method and for each mock BBH in the sample,its mass,mass ratio and GW frequency at the starting time (fi) of DO observations can be obtained.Considering that DO is extremely sensitive for detecting BBHs,we first consider a half year of observations and generate a mock BBH sample that can be detected by DO with an S/N (ϱ) threshold of 8.

The S/N of each BBH is given by Maggiore (2008)

with

Note here θ,φ,ψ and ι are the polar angle,azimuthal angle,polarization angle and inclination angle of the source in the detector frame respectively.After averaging over all possible directions and inclinations,we have

Combining Equations(20),(21)and(25),the averaged S/N is

where j refers to detector“1”(with α0=0,Φ0=0)or“2”(with α0=π/4,Φ0=0).Given Γab,then we have

and thus we can estimate the uncertainties in the measurements of Ξ as

In particular,the angular resolution ΔΩ is de fined as

For both DO detectors,the GW strain signal can be described by Equation (21).The GW strain phase evolution includes the polarization modulation and Doppler modulation

The motion of DO detectors around the helio-center causes Doppler modulation of the GW phase as

where R=1 au is the orbital radius of the detector and

With the above formalism,we can estimate S/N values for all the mock BBHs that can be detected(ϱ ＞8)by DO within a half year observation and the precision of parameter measurements for each “detectable” mock BBH.For all those“detectable” BBHs,Figure 8 plots the distributions of their S/N (bottom-right panel),redshift (top-right panel),measurement errors of luminosity distance (top-left panel) and localization (bottom-left panel).Apparently more than hundreds of mock BBHs can be detected with ϱ ≿30,which basically have accurate determinations of their luminosity distances and locations.The dashed histograms in this figure show the distributions of those mock BBHs (totally 90 and typically at z ≾2) that haveσdL＜0.02 and ΔΩ≾0.1 deg2.This subsample may be used to do cosmological inference as we have demonstrated in Section 4 that the mock BBHs with precise luminosity distance measurements and localization can be used to obtain tight constraints on the Hubble constant and other cosmological parameters.

Figure 8.Distributions of mock BBH mergers expected to be detected by DO.Solid and dashed histograms display the distribution of the whole mock sample of BBHs “detected” by DO with S/N ＞8 in a half year observation period(totally 10 269)and that of a subsample with ＜0.02(totally 90)adopted for cosmological parameter inference,respectively.

Figure 9.Posterior probability distributions of h and Ωm obtained from the mock BBH mergers“detected”by DO with luminosity distance uncertainties of δdL/dL＜0.02.

Note here that we do not consider the bias correction introduced in Section 4.2 as the DO sample we adopted above all have precise measurements on the luminosity distances(δdL/dL＜0.02) and localization areas (δΩ ≾0.01 deg2).If including sources with larger luminosity distance errors and localization errors,one may need to consider the bias correction.

For comparison,the accuracy of recent h measurements by using the local distance ladder is about 2% level (e.g.,Riess et al.2019).Surveys such as those with Euclid,the Square Kilometre Array (SKA),eBOSS and DESI are expected to reach sub-percent precision on Baryon Acoustic Oscillation(BAO) measurements of the Hubble parameter at redshift 0.1＜z＜3 and give constraints on h close to 1% in the best cases (Wang et al.2017; Bengaly et al.2020).As for GW standard sirens,Chen et al.(2018) predict a 2% measurement of h with Advanced LIGO-Virgo within five years,which is expected reach about 0.3% for future ET with 10 years of observation (Zhang et al.2019).In addition,with five-years of observation by future space-based GW detector Taiji,the constraint on h is expected to be～1.3%using massive BBH as standard sirens and～1.0%with the network of LISA and Taiji(Wang et al.2021).This suggests that,with the observation of DO,stellar mass BBHs as dark sirens,as an independent measurement of cosmological parameters,is quite promising.

6.Conclusions

GW events can be used as standard sirens to constrain cosmological parameters if their redshift can be measured by EM observations.Although a large fraction of GW events may not have EM counterparts and cannot have direct redshift measurements,e.g.,BBH mergers,one may still use it as dark standard sirens to statistically constrain cosmological parameters provided that the redshifts of their host candidates can be obtained from deep galaxy surveys.In this paper,we investigate this dark siren method in detail by using mock GW events and galaxy catalogs.We find that the Hubble constant can be constrained well(on percentage level or better)by using a few tens or more of BBH mergers at redshift up to 1 with accurate luminosity distance measurements (relative error≾1%) and localization (sky coverage ≾0.1 deg2),though the constraint may be significantly biased if these measurements are less accurate and the bias cannot be removed by simply increasing the number of adopted GW events.We analyze this systematic bias in detail and introduce a correction method to remove it.We demonstrate that this correction method works well in most cases,though it may not work well when the measurement errors of GW luminosity distances and localization areas are too large.Our results suggest that GW events as standard dark sirens can be utilized to robustly constrain the cosmological parameters once their luminosity distances and localization areas are well determined.So,selection of sources with sufficiently high precision measurements of the luminosity distances and localization areas is important for obtaining robust constraints on the cosmological parameters.

Finally,we simulate future observations of BBH inspirals and mergers by DO,according to current constraints on BBH merger rate and distributions of BBH properties.We find that DO in a half year observation period may detect about one hundred BBHs with S/N ϱ ≿30,relative luminosity distance error ≾0.02 and localization error＜0.01 deg2.We predict that the Hubble constant can be constrained to the～0.1%–1%level using these DO BBHs as dark sirens,which is comparable with the accuracy of other methods in the near future.This dark siren method,independent of those methods that rely on the local distance ladder,BAO and CMB observations,would provide insight into the tension between the early-and late-Universe measurements of the Hubble constant.We also demonstrate that the constraint on the Hubble constant using this dark siren method is robust and does not depend on the choice of the prior for the properties of BBH host galaxies.

Acknowledgments

This work is partly supported by the National Key Program for Science and Technology Research and Development(Grant Nos.2020YFC2201400,2020SKA0120102 and 2016YFA0400704),the National Natural Science Foundation of China (Grant No.11690024),and the Strategic Priority Program of the Chinese Academy of Sciences (Grant XDB 23040100).