Constraining the cosmological parameters using gravitational wave observations of massive black hole binaries and statistical redshift information

Space-borne gravitational wave detectors like TianQin are expected to detect GW signals emitted by the mergers of massive black hole binaries. Luminosity distance information can be obtained from GW observations, and one can perform cosmological inference if redshift information can also be extracted, which would be straightforward if an electromagnetic counterpart exists. In this paper, we concentrate on the conservative scenario where the EM counterparts are not available, and comprehensively study if cosmological parameters can be inferred through a statistical approach, utilizing the non-uniform distribution of galaxies as well as the black hole mass-host galaxy bulge luminosity relationship. By adopting different massive black hole binary merger models, and assuming different detector configurations, we conclude that the statistical inference of cosmological parameters is indeed possible. TianQin is expected to constrain the Hubble constant to a relative error of about 4%-7%, depending on the underlying model. The multidetector network of TianQin and LISA can significantly improve the precision of cosmological parameters. In the most favorable model, it is possible to achieve a level of 1.7% with a network of TianQin and LISA. We find that without EM counterparts, constraints on all other parameters need a larger number of events or more precise sky localization of GW sources, which can be achieved by the multidetector network or under a favorable model for massive black hole mergers. However, in the optimistic case, where EM counterparts are available, one can obtain useful constraints on all cosmological parameters in the Lambda-CDM cosmology, regardless of the population model. Moreover, we can also constrain the equation of state of the dark energy without the EM counterparts, and it is even possible to study the evolution of EoS of the DE when the EM counterparts are observed.

Space-borne gravitational wave detectors like TianQin are expected to detect gravitational wave signals emitted by the mergers of massive black hole binaries.Luminosity distance information can be obtained from gravitational wave observations, and one can perform cosmological inference if redshift information can also be extracted, which would be straightforward if an electromagnetic counterpart exists.In this paper, we concentrate on the conservative scenario where the electromagnetic counterparts are not available, and comprehensively study if cosmological parameters can be inferred through a statistical approach, utilizing the non-uniform distribution of galaxies as well as the black hole mass-host galaxy bulge luminosity relationship.By adopting different massive black hole binary merger models, and assuming different detector configurations, we conclude that the statistical inference of cosmological parameters is indeed possible.TianQin is expected to constrain the Hubble constant to a relative error of about 4%-7%, depending on the underlying model.The multidetector network of TianQin and LISA can significantly improve the precision of cosmological parameters.In the most favorable model, it is possible to achieve a level of 1.7% with a network of TianQin and LISA.We find that without electromagnetic counterparts, constraints on all other parameters need a larger number of events or more precise sky localization of gravitational wave sources, which can be achieved by the multidetector network or under a favorable model for massive black hole mergers.However, in the optimistic case, where electromagnetic counterparts are available, one can obtain useful constraints on all cosmological parameters in the Lambda cold dark matter cosmology, regardless of the population model.Moreover, we can also constrain the equation of state of the dark energy without the electromagnetic counterparts, and it is even possible to study the evolution of equation of state of the dark energy when the electromagnetic counterparts are observed.
The key to the success of the standard siren method involves obtaining redshift information for GW sources.One needs to either (1) identify the EM counterpart and the host galaxy of GW sources, or (2) obtain the statistical redshift distribution of candidate host galaxies, based on our knowledge of galaxy clustering [20].Novel methods have also been developed to obtain the redshift information of GW sources, using such effects/information as [30]: (3) the strong gravitational lensing of GWs [31][32][33]; (4) the mass distribution function of compact binary GW sources, such as [34,35] for BNS and [36,37] for stellar-mass binary black hole (StBBH), dependent on our understanding of the history of binary mergers; (5) the redshift probability distribution of compact binary mergers based on their intrinsic merger rates [38], requiring a large number of GW events and also dependent on our understanding of the history of binary mergers; (6) the phase correction of BNS merger GWs due to tidal effects [39][40][41], requiring detectors with higher sensitivity such as third-generation ground-based GW detectors; (7) the evolution of GW phase with cosmological expansion [42][43][44], requiring detectors with both higher sensitivity and operating in the decihertz band.
Third-generation ground-based GW detectors such as the Einstein Telescope (ET) [54] and Cosmic Explorer (CE) [55] are expected to detect thousands of GW events with a redshift concentration at z ≈ 2 and a horizon of z ≈ 10 [35,38,56,57].The high-redshift GW events detected by ET or CE can provide precise measurements of H 0 and other cosmological parameters, such as the density of dark matter Ω M , and of dark energy Ω Λ [35, 37-39, 52, 58-62].In addition, a large number of high-redshift GW events would allow us to reconstruct the dark energy equation of state (EoS) and expansion dynamics by nonparametric methods [59,63,64].
Space-borne GW detectors like TianQin will open for exploration the millihertz to Hertz band of the GW spectrum.TianQin is a constellation of three satellites orbiting around the Earth, using drag-free control to lower noise and measure GW effects through laser interferometry [65,66].TianQin is expected to observe multiple different types of GW source, including StBBH inspirals [67], extreme mass ratio inspirals (EMRIs) [68], Galactic compact binaries [69] and massive black hole binary (MBHB) mergers [70].Such detections are also expected to put stringent constraints on deviations from general relativity or testing specific gravity theories [71,72].Studies have revealed the availability of stable orbits that fulfils the requirement for GW detections [73][74][75][76].Moreover, the technology demonstration satellite, TianQin-1, has met the design requirements [77].
In this paper, we focus on the ability of space-borne GW detectors to constrain the cosmological parameters using MBHB merger GW signals.Some studies have suggested that MBHB mergers may have observable electromagnetic signatures [95][96][97][98], and some thus assume the availability of an EM counterpart when performing GW cosmology studies [78,81,83,85,99,100].In order to be conservative, however, we set as our default assumption that no EM counterparts are detectable for our GW sources, so that we rely on the luminosity distance information from the GW detections combined with statistical information about their host galaxy redshifts to constrain the cosmological parameters.Our method can be simply described as follows: GW detection provides a localisation error cone for the three-dimensional (3D) position of the GW source; one can use the redshift distribution of the galaxies within the cone as the proxy for the GW redshift.However, there could be a lot of galaxies within the localisation error volume that will cause contamination and limit the effectiveness of the dark standard siren method.We therefore study how using the relationship between the central massive black holes (MBHs) and their host galaxies [101][102][103][104] can better pinpoint the host and improve the cosmological Throughout this paper, we consider a Universe in which the spacetime can be described by the Friedmann-Lemaître-Robertson-Walker metric, and adopt the Lambda cold dark matter (ΛCDM) model as our fiducial model, with the dark energy EoS being described by a constant ω ≡ p Λ /ρ Λ = −1.In this model the expansion of the Universe can be characterized by the Hubble parameter H(z) ≡ ȧ/a, with a being the scale factor, the current value of which is a 0 , and with redshift z defined as z + 1 ≡ a 0 /a.The Hubble parameter can therefore be expressed as where the Hubble constant H 0 ≡ H(z = 0) describes the current expansion rate of the Universe, and Ω M , Ω K and Ω Λ are respectively the current dimensionless fractional densities for the total matter, curvature and dark energy with respect to the critical density.They satisfy the relationship Ω M + Ω K + Ω Λ ≡ 1.
We will also consider the possibility that the dark energy component has dynamical properties, by adopting the Chevallier-Polarski-Linder (CPL) parametrization [106,107] This yields the equation Estimates of the cosmological parameters can be inferred by fitting the relationship between observed luminosity distances D L and redshifts z.This relationship encodes information about the expansion history of universe, and is given by where c is speed of light in vacuum.We adopt cosmological parameters derived from various recent observations, like Planck [24], H 0 = 67.8km/s/Mpc, Ω M = 0.307, Ω Λ = 0.693 and we adopt w 0 = −1, w a = 0, respectively, consistent with the observed galaxy distribution [108].

B. Standard sirens
For the inspiral of a compact binary system with component masses m 1 and m 2 , the frequency domain GW waveform can be expressed as [109] where G is the gravitational constant, M = η 3/5 M is the chirp mass, η = m 1 m 2 /M 2 is the symmetric mass ratio, M = m 1 + m 2 is the total mass, and Φ(f ; M z , η) is phase of the waveform depending on the parameters M and η.The chirp mass M largely determines the overall evolution of the GW waveform, but the parameter directly measured from the data is actually the redshifted chirp mass M z ≡ M(1 + z).
One can see from Eq. ( 5) that the luminosity distance D L of the binary has a direct impact on the measured waveform amplitude.Therefore, GW observations of compact binary coalescences can be used to estimate their corresponding luminosity distance directly.If the redshifts of such mergers can be inferred through other observational or theoretical channels, one can use information from both luminosity distance and redshift to constrain the cosmological parameters [20].Such a one-stop measurement of luminosity distance makes the inspiral signal of compact binary systems desirable objects for cosmological studies, since they are largely immune to systematic errors caused by intermediate calibration stages like those required by type Ia supernovae; consequently, they have been coined as "standard sirens" .
There is a catch, however.Both D L and z are needed in order to use Eq. ( 4) to infer cosmological parameters.Moreover, as noted above, the redshift z is deeply intertwined with the chirp mass M and can not be solely determined by GW observations of the inspiral.The measurement of redshift information thus relies on extra information, like the identification of the host galaxy or the direct observation of an EM counterpart.

C. Bayesian framework
We adopt a Bayesian framework to infer cosmological parameters through the GW observations of massive black hole binary mergers.Consider a set of data composed of N GW observations, D ≡ {D 1 , D 2 , ..., D i , ..., D N }, as well as an EM data set S derived from EM observations.Then one can derive the posterior probability distribution of the cosmological parameters Ω as where I indicates all relevant background information.Since the normalisation factor, also known as the Bayesian evidence, p(D, S|I), is irrelevant in the calculation, the posterior can be written as We can classify parameters into three categories: the common parameters (which affect both GW waveforms and EM observations); GW-only parameters θ , and EM-only parameters φ .Throughout this paper, we identify the common parameters as the luminosity distance D L , the redshift z, the longitude α, the latitude δ, the total mass M of the compact binary, and the bulge luminosity L bulge of the host galaxy.We can express the likelihood as where M z = M (1 + z) is the redshifted total mass of the GW source.Notice that we introduce a correction term β( Ω|I) to eliminate the effect of selection biases [14,49,110].The integrand in the numerator of Eq. ( 8) can be factorized as For details on the derivation of Eq. ( 9), please refer to Appendix A. The general mathematical treatment of the likelihood p(D i |D L , α, δ, M z , θ , I) can be given by [111] p where •|• is the inner product as defined in Eq. (22).When the host galaxy of the GW signal cannot be uniquely identified, one can set p(S|z, α, δ, L bulge , φ , I) as a constant [49].We assume that p 0 (D L |z, Ω, I) ≡ δ D L − DL (z, Ω) depends on the cosmological model and p 0 (M z |z, L bulge , Ω, is based on the relation of massive black hole mass with galactic luminosity [102,103], where the total mass M (L bulge ) is a function of the galactic bulge luminosity L bulge .
The EM measurements of parameters like position of the galaxies (α, δ) are much more precise compared with the measurements from GW observations.The prior in Eq. ( 9) can therefore be approximated as [49,51,87] where N gal is the number of galaxies in our EM catalog, and N (x|x, σ x ) is a Gaussian distribution on x, with expectation x and standard deviation σ x , here x = {z, L bulge }.However, if the luminosity function Φ(L) is also regarded as redshift-dependent (as in [112][113][114] etc.), then Eq. ( 11) needs to be supplemented, as shown in Eq. (31).
Marginalizing over the parameters D L , M z , θ , and φ , Eq. (8) becomes We next examine the correction term β( Ω|I).Bias of the inferred cosmological parameters might arise from two sources: the GW data and the EM information.The EM information contains the following bias: (1) the 3D error volume is cone like, so this will in general lead to a bias towards higher redshift due to the larger associated volume; and (2) the incompleteness of the catalog will introduce a Malmquist bias, where brighter galaxies are disproportionately recorded and weighted [115].To obtain cosmological constraints from "dark" standard sirens depends on the non-uniform distribution of galaxies, due to large scale structures (LSS) or smaller-scale clustering of galaxies, so the correction term should exclude the influence of LSS and galaxy clustering information as far as possible.We account for the aforementioned two biases by counting the detectable galaxies over the whole sky, and considering the redshift evolution of this galaxy count.We assume that the EM observations are isotropic within the survey region, and define the prior for the redshift distribution of galaxy catalog as where ∆z is chosen to be much larger than the typical redshift scale of the LSS.The correction term after marginalization is then Notice that the GW selection effect is accounted for by integrating over only events detectable by GW detectors.
In addition, if the catalog of survey galaxies contains two or more sky areas with different observation depths, then p c (z| Ω, I) and β( Ω|I) terms need to be calculated separately for each sky area.For our analysis, the incompleteness of a galaxy catalog can introduce two effects.First, it could bring Malmquist bias so that more distant galaxies are disproportionately weighted.Such a systematic bias is caused by the selection effect that a catalog tends to be more complete for brighter galaxies.The correction term β( Ω|I) accounts for this bias.On the other hand, a less complete catalog means there is a higher chance for the host galaxy cluster to be missed.In this case, an additional bias could be introduced to the analysis.In Sec.V E we demonstrate that such a large deviation can be identified through application of a consistency check.Meanwhile, in [53] the authors take a step further than Eq. ( 14), so that more distant events are down-weighted and a consistency check is considered unnecessary.

D. Parameter estimation of the GW sources
We define the sensitivity curve S n (f ) in terms the expected power spectral density S N (f ) according to the relation, where T (f ) is the sky and polarization averaged response function of the detector, and Here L is the arm length, S a is the acceleration noise, S x is the positional noise, and S c (f ) is the Galactic foreground noise.For TianQin, S .In addition, we added the Galactic foreground on top of the sensitivity curve assuming an observation time of 4 years according to [105].On the other hand, no foreground is added for TianQin and it has been suggested that the anticipated foreground will be below TianQin's sensitivity [69,116].The sky-and polarization-averaged sensitivity curves of TianQin and LISA are shown in Fig. 1.A simultaneous observation from multiple detectors can improve the sky localisation of the GW source.Specifically, the long baseline between different GW detectors makes it possible to use the time delay information to perform sky localisation [5-7, 14, 86, 117].
Suppose that the kth detector record the GW strain h k (t), the strain takes the form where ι is the source inclination angle, A(t) and Φ(t) are the amplitude and phase of the GW respectively, F k +,× (t) are the response functions of the kth detector, Φ k D (t) is the Doppler frequency modulation due to the Doppler effect of the solar orbital revolution of the detector, and Φ k P (t) is the phase modulation.The expressions for the phase modulation and the response function are closely related to the orbit of the detectors.We follow [70,118] for TianQin, and [119,120] for LISA, adopting the low-frequency limit.
For a circularised binary black hole merger event, the GW signal can be generally described by a set of parameters including the binary component masses m 1 and m 2 (which can be reexpressed as chirp mass M and symmetric mass ratio η), the luminosity distance D L , the longitude α and latitude δ of the source sky position, the coalescence phase φ c , the coalescence time t c , and the direction of the orbital angular momentum of the GW source (α L , δ L ) (which can be reexpressed as an inclination angle ι and polarization angle ψ).More parameters would be needed if the spins of the black holes are included.
Let us consider a detector network including N independent detectors, for which the frequency domain GW signal h(f ) can be written as where h k (f ) denotes the Fourier transform of h k (t).The signal-to-noise ratio (SNR) ρ of a GW signal h(f ) can be defined as [111,121] where the inner product symbol •|• is defined as where * represents complex conjugate, Re denotes the real component, and S k n (f ) are the sensitivity curve functions of the kth detector respectively when the response functions F +,× (t) adopt the low-frequency approximation [70,105].
For a GW signal from a binary characterised by physical parameters θ = (M, η, D L , α, δ, cos ι, ψ, φ c , t c ), the inverse of the corresponding Fisher information matrix (FIM) sets the Cramér-Rao lower bound for the covariance matrix [122].The FIM can be written as where θ m indicates the mth parameter.With Σ = Γ −1 , we take the estimation uncertainty on parameter θ m as ∆θ m = √ Σ mm , with the sky localization error ∆Ω given by the combination ∆Ω = 2π| sin δ| Σ αα Σ δδ − Σ 2 αδ .If the GW signals are gravitationally lensed, or the peculiar velocity of its host galaxy is not properly accounted for, the inferred luminosity distance would be in error.Therefore, in addition to the measurement error from GW observation, weak lensing and peculiar velocity can both contribute to the intrinsic uncertainty of the estimated D L .We denote σ GW D L as the uncertainty arising from the GW observation, and σ tot D L as the overall uncertainty including all sources of error.We adopt a fitting formula from [123] (also see [124]) to estimate the weak lensing error σ lens D L , given by where C l = 0.066, β l = 0.25 and α l = 1.8.For the error, σ pv D L , due to peculiar velocity we adopt the fitting formula from [125,126], with v 2 = 500 km/s as the root mean square peculiar velocity of the host galaxy with respect to the Hubble flow [127].We then define the total uncertainty Throughout this paper, we use σ tot D L for our likelihood calculation.However, we remark that the above is a conservative estimate of the total uncertainty, since "delensing" methods like the use of weak lensing maps [128], observations of the foreground galaxies [129] or deep shear surveys [130] could in principle be used to alleviate σ lens D L , while peculiar velocity maps could also be used to reduce σ pv D L [14,131].We do not apply such corrections in this paper; hence we can expect better results for GW cosmological inference from realistic future detections.

III. SIMULATIONS
For the purpose of our paper, we require a catalog of massive black hole mergers that mimic our understanding of the real Universe.In this section, we describe the simulation of these catalogs, and how we use the 3D localisation information derived from a GW detection, as well as the empirical M MBH − L bulge relation, to allocate probabilities to candidate host galaxies.

A. Massive black hole binary mergers and galaxy catalog
Following previous studies, e.g., [70,132], we adopt the massive black hole binary merger populations from [133].Both the "light-seed" scenario and the "heavy-seed" scenario are considered for the seeding models of massive black holes.In the light-seed scenario, seed black holes are assumed to be the remnants of first generation (or population III) stars and the mass of the seed black hole is around 100 M [134,135].This model is later referred to as popIII.In the heavy-seed scenario, the MBH seeds are assumed to be born from the direct collapse of protogalactic disks that may be driven by bar instabilities [132], with a mass of ∼ 10 5 M .The critical Toomre parameter Q that determines when the protogalactic disks become unstable is set to 3 [136].Two models are derived from this scenario, with Q3d considering the time lag between the merger of MBHs and merger of galaxies, and Q3nod, which does not consider such a time lag.In both scenarios, the seed BHs grow via accretion and mergers, eventually becoming the massive black holes, while the evolution of MBHs are deeply coupled to the evolution of their host galaxies [137].
In addition to the simulated catalog of GW mergers, we also require a simulated galaxy catalog, so that we can associate each MBHB merger with a certain host galaxy and use multimessenger information to infer cosmological parameters.TianQin has the ability of observing very distant mergers, but no existing galaxy survey project can extend to these distances so we choose to adopt the MultiDark Planck (MDPL) cosmological simulation [108] from the Theoretical Astrophysical Observatory (TAO) [138] for this purpose.The MDPL simulated a catalog of galaxies based on an N -body simulation which tracks the evolution of dark matter halos assuming a Planck cosmology [139], with 3840 3 particles and a box side length to h −1 Gpc (where h ≡ H0 100 km/s/Mpc ), and the simulation is performed from z = 100 to z = 0. Additional information such as galaxy evolution was obtained from the semi-analytic galaxy evolution model [140], and the luminosity distribution from Conroy et al. [141].

B. GW event catalog
We first simulated MBHB mergers according to the three models (popIII, Q3d, and Q3nod), where only events with network SNR ρ > 8 are later used.Notice that although TianQin has the ability to observe very distant events, it is anticipated that complete galaxy catalogs would be extremely hard to obtain for galaxies beyond redshift z = 3, therefore we apply a redshift cut beyond this point.The remaining parameters are then generated from simple distributions: all angle parameters are chosen uniformly in solid angle, , and we choose the merger time t c ∈ U[0, 5] years, merger phase φ c ∈ U[0, 2π] and the spin χ 1,2 ∈ U[−1, 1].Based on these parameters, we generate the GW waveform from the IMRPhenomPv2 model [142].
Throughout this analysis, we consider multiple configurations for the space-borne GW detectors as described below: • (i) TianQin: the default case where three satellites form a constellation and operate in a "3 month on + 3 month off" mode, with a mission life time of 5 years [65]; • (ii) TianQin I+II : a twin constellation of satellites are used that have perpendicular orbital planes to avoid the 3-month gaps in data [67,70]; • (iii) LISA: we consider the mission as described in [80,105], with a nominal mission life time of 4 years; • (iv) TianQin+LISA: TianQin and LISA observing together, with 4 years of overlap in operation time; • (v) TianQin I+II+LISA: similar to above but with the TianQin I+II configuration considered.
Previous work indicates that both TianQin [70] and LISA [132] have good detection capabilities for MBHB mergers.In Table I we list the anticipated total detection rates for mergers at redshift z < 3 under different detector configurations.It is worth noting that the actual merger rate of MBHB is likely to increase by about twice as much as the three population models predict [132], and our actual detection rate would therefore also likely increase by about twice as much.In order to identify the host galaxy, the most important information that we gather from the GW observation would be the spatial localisation, equivalently α, δ, and D L .In Fig. 2, we illustrate the marginalised distribution on the sky localisation error ∆Ω, as well as the relative error on the luminosity distance σ D L /D L .
As shown in Fig. 2, a typical MBHB merger can be localised to better than 1 deg 2 with TianQin alone, while a combination of both TianQin and LISA can improve the localisation precision by a factor of 1 ∼ 2 orders of magnitude, where a small fraction of sources can even be localised to within 10 −4 deg 2 .However, this is generally not sufficient to pinpoint the host galaxy uniquely, especially considering the relatively large uncertainty in luminosity distance as well as the large distance range to which the detectors can reach.
Thanks to the relatively high SNR, the typical D L uncertainty arising from the GW observations is around the level of 1% across the three models, while the combination of both TianQin and LISA can further improve the precision by a factor of about three.However, as indicated in the bottom left panel of Fig. 2, often the GW measurement error is the subdominant contribution.Therefore the overall error in luminosity distance, after considering the effects of weak lensing and peculiar velocity, is of order 10%.

C. MMBH − L bulge relation
Once the MBHB mergers and galaxy catalogs are ready, the next step is to associate the appropriate galaxy as the host galaxy for each merger.Throughout this paper, we link the MBHB mergers and their host galaxies through the M MBH − L bulge relation.We adopt the form of this relation as [101][102][103] where M 0 , k c and L 0 are fitting parameters from astronomical observations.In the K-band, the values of the fitting parameters are log 10  27) to describe the M MBH − L bulge relation, it is accompanied with significantly larger scatters [104].Therefore, in the low-mass end, we adopt a different intrinsic scatter σ int M −L ≈ 0.5 [104].This means that the intrinsic scatter for MBHB mass from the M MBH − L bulge relation is a step function And in the following, we will refer to the luminosity in terms of solar luminosity and black hole mass in terms of solar masses, denoted by L and M , respectively.

D. Simulation of the statistical redshift information
The dark standard siren method relies on the measurement of luminosity distance from GW observations, and the inference of redshift information from galaxy catalogs.So the appropriate identification, or at least the association, of the host galaxy with the merging MBHB is of the utmost importance.
We first obtain a conservative range of luminosity distance for the possible host galaxy, where DL is mean estimated value.Notice that we shall not use the galaxy luminosity distance information directly, otherwise it is pointless to introduce GW observations for constraining cosmology.
This choice of parameters ensures that mainstream estimates of Hubble constant, albeit in disagreement with each other, are nonetheless well encapsulated by our prior [22,24,[26][27][28].
It is worth noting that the spatial localisation errors of the GW sources are indeed of irregular shape, instead of 3D ellipsoids or cylinders as some previous analyses have assumed for simplicity.We illustrate the localisation error for a typical event in Fig. 3. Notice also that in order to further simplify the ensuing calculation, we replace the elliptic sky localisation error with a circular shape but of the same area, which shall not alter the statistical conclusions of our paper.
We then artificially assign a random galaxy to be the actual host galaxy.For a galaxy with bulge luminosity L bulge , we assign a weight according to a log-normal distribution with the M MBH − L bulge relation.Next, we aim to simulate the real observation, to obtain a catalog of galaxies and to assign a probability of each galaxy hosting the merging black hole binary.Furthermore, we simulate the Malmquist bias to be more realistic.For the jth galaxy with luminosity L j and redshift z j , the probability of it being recorded is where log 10 L limit (z j , Ω) = (+4.83−mlimit −5)

2.5
+ 2 log 10 ( ), dependent on a given set of cosmological parameters Ω, and we adopt a limiting magnitude of m limit = +24 mag and equal measurement error of luminosity σ log 10 L = 0.04 (correspond to an uncertainty in magnitude of 0.1 mag) [143,145,146,171].
In the analysis stage, we try to eliminate the Malmquist bias by introducing a redshift-dependent correction factor using the luminosity function of both the galaxies [112,113,147,148] and the bulge.We need to manually augment sample sizes for further galaxies.We first divide the error box into multiple small regions both in sky location as well as in redshift.For each region, we can calculate the number of supplementary galaxies Nsup using a luminosity function Φ(L), Nsup = N obs where we derive the luminosity function Φ(L) from the MDPL simulated catalog [141].For a given GW source, the debiased prior of the location is determined by the possible host galaxies, which can be expressed as where N tot = (N obs + Nsup ), and zobs , ᾱobs , and δobs are the mean values of redshift, longitude and latitude for the observed galaxies in the small region.For nearby galaxies (z < 1), we assume that spectroscopic redshift information is available and therefore the error can be neglected [149].However, for further galaxies the redshift if more likely obtained through photometric measurement, which is assumed to be associated with an error ∆z = 0.03(1 + z) [150,151].Therefore, we adopt The bulge luminosity function with an apostrophe, Φ (L bulge ), is defined as where L min bulge is the minimum bulge luminosity of observed galaxies in the small region, and we derive the bulge luminosity function Φ(L bulge ) from the MDPL simulated catalog.We remark that the second term in Eq. (31) does not represent a new batch of galaxies, but rather an adjustment to the weights of existing galaxies.
Next we assign different weights for galaxies with different positions and bulge luminosities.We consider the two following methods: • (i) fiducial method : each galaxy in the spatial localisation error box of the GW source has equal weight regardless of its position and luminosity information; • (ii) weighted method : the weight of a galaxy is the product of both its positional weight and bulge luminosity weight.
The positional weight is simply determined by the 3D space localisation from the GW parameter estimation.The luminosity-related weight, on the other hand, is more complicated.For an observed galaxy, the weight is assigned through a log-normal distribution, with the expected redshifted mass value being the mean value derived from the GW parameter estimation, and standard deviation σ log 10 M .And for a manually supplemented galaxy, a further integration of this log-normal distribution over luminosity up to L min bulge is needed.A detailed expression for both the positional weight and the luminosity-related weight for a given galaxy is provided in Appendix B.
The luminosity-related weight of galaxies can be affected by the intrinsic scatters of the M MBH −L bulge relation as well as the uncertainties related to the measurements.We define σ log 10 M as the root sum squared of the intrinsic scatter of the relation between the central MBH mass and the galactic bulge luminosity σ int M −L , the measurement error on the total mass of GW sources σ GW log 10 M , and the error related to the bulge luminosity k c σ EM log 10 L bulge .We adopt σ GW log 10 M ≡ 0.05 as a conservative choice since the mass parameter can usually be accurately determined through GW observation.For the measurement error on the bulge luminosity, since we do not have complete information on the bulge luminosity distribution for high redshift galaxies, we adopt the galaxy luminosity error as a proxy through σ EM log 10 L bulge We demonstrate the effect of the debias and different weighting in Fig. 4. For the top panel, we do not correct for Malmquist bias, and there is an apparent excess of galaxies at low redshift.In the middle and bottom panel, the Malmquist bias is corrected by manually introducing supplementary galaxies; in these cases the distribution roughly follows the prior.For the top and middle panel, all galaxies within the error box are assumed equally likely to be the host galaxy of the MBHB merger, while for the bottom panel, we assign weights to the galaxies using the weighted method.One can observe that the bottom panel is less smooth, and the injected source stands out from the contaminating galaxies.

IV. COSMOLOGICAL CONSTRAINTS
In this section, we explore quantitatively the prospects for parameter estimation with GW cosmology.We perform a series of Markov chain Monte Carlo (MCMC) studies using the widely adopted library emcee, a Python package that implements an affine-invariant MCMC ensemble sampler [152,153].For the detector(s), we consider various scenarios, including: TianQin, TianQin I+II, LISA, TianQin+LISA, and TianQin I+II+LISA.We investigate both the dark standard siren case, where no EM counterpart is available, and the bright standard siren case where we can use an EM counterpart to identify the host galaxy and therefore to obtain explicitly the redshift from EM observations.For the popIII and Q3d models we generate 200 mock GW events from MBHB mergers, while for the Q3nod model, the event rate of which is expected to be higher than the other two models, we generate 600 mock events.These events form three sample pools, and the GW events on which each realization of the cosmological parameter estimation relies are chosen randomly from these three pools.
In order to comprehensively probe the systematic and random errors in the inferred cosmological parameters, we repeat the cosmological analysis multiple times with independent runs for each detector configuration/MBHB population model/weighting method.

A. TianQin and TianQin I+II
We first consider the most pessimistic scenario, where TianQin is operating alone and no EM counterpart is expected to be observed.For convenience, in this analysis we simply select all galaxies that fall into the range of 9∆Ω × [z min , z max ] (although the actual error box would correspond to a shape type characterised by Fig. 3).
As illustrated by [154], the SNR accumulates rather quickly just before the final merger.Therefore, the difference in the precision of cosmological parameter estimation between TianQin and TianQin I+II is rooted in their different detection numbers.The distributions of the constraint precisions of various cosmological parameters for TianQin (left column) and TianQin I+II (right column) under the three MBHB population models are shown in Fig. 5.In order to eliminate the random fluctuations caused by the specific choice of any event, we repeat this process 24 times, and the constraint results are presented in the form of boxplots.
The top, middle, and bottom rows represent the estimation precision of H 0 , Ω M , and Ω Λ , respectively.The constraining ability of these parameters are decreasing as with this order.In each panel, from left to right, we show the results adopting popIII, Q3d, and Q3nod as the underlying astrophysical model, respectively.In each model, we generate a number of GW event catalogs, assuming an event rate that follows a Poisson distribution with the rate parameter determined by Table I.It can be observed that more events leads to better constraints on the cosmological parameters.Furthermore, the weighted method (red box) shows better constraining ability than the fiducial method (gray box).
In general, under the models with a lower MBHB event rate, such as popIII and Q3d models, we can only constrain the Hubble constant, which is almost entirely determined by low redshift GW events.The constraints of other cosmological parameters require more GW events, which is feasible under the Q3nod model.For TianQin, the relative precision on H 0 is estimated to be 7.8%, 7.5% and 4.2% via the fiducial method; by using the weighted method, the relative precision can reach 6.9%, 6.5% and 3.3%, for the popIII, Q3d and Q3nod models, respectively.For TianQin I+II, the number of MBHB mergers would be boosted by a factor of about 2, and the relative error of H 0 reduces to 6.0%, 6.0%, and 2.0% using the weighted method, respectively.Under the Q3nod model, Ω M and Ω Λ are expected to be constrained to a relative precision of 34.9% (25.1%) and 26.1% (21.5%) for TianQin (TianQin I+II), respectively.
Finally, we study the constraining power of the GW observations on the EoS of dark energy, assuming the CPL model for its redshift evolution.Here we fix H 0 = 67.8km/s/Mpc, Ω M = 0.307, and Ω Λ = 0.693.We find that, for all astrophysical models and detector configurations, w a can hardly be constrained; however, meaningful constraints can be obtained on w 0 , as shown in Fig. 6.And similarly, the weighted method also leads to more precise constraints on the parameters of EoS of dark energy compared to the fiducial method.If using the weighted method, the relative precision of w 0 can be reach a level of 36.7%, 37.2%, and 13.8% for TianQin, and 26.6%, 29.6%, and 8.1% for TianQin I+II -for popIII, Q3d, and Q3nod models respectively.

SRS,,,4G4QRG
7LDQ4LQ,,, SRS,,,4G4QRG SRS,,,4G4QRG FIG. 5. Boxplot of the precision of the estimated cosmological parameters for the three MBHB models, assuming two detector configurations, i.e., TianQin (left column) and TianQin I+II (right column).The top, middle, and bottom rows illustrate results for three cosmological parameters, H0, ΩM , and ΩΛ, respectively.The horizontal gray dashed line represents a fiducial 68.27% statistical interval from the prior.For each result, the box represents the range of 25% − 75% of the data distribution, the upper limit of the whisker length is 1.5 times the box length, and the crosses are outliers.
In each box, the dot represents the mean value, the short horizontal line represents the median value, using the fiducial method (gray) and the weighted method (red), respectively.

SRS,,,4G4QRG
6 × 10 1 4 × 10 1 2 × 10 1 FIG. 6. Boxplot of the precision of the estimated dark energy EoS parameters for the three MBHB models, assuming two detector configurations, i.e., TianQin (left column) and TianQin I+II (right column), and keeping non-CPL parameters fixed.The top and bottom illustrate results for three cosmological parameters, w0 and wa, respectively.The horizontal gray dashed line represents a fiducial 68.27% statistical interval from the prior.The gray and red boxplots represent the distribution of estimation precisions of the parameters using the fiducial method and the weighted method, respectively.
Besides, in order to more graphically represent the advantages of the weighted method over the fiducial method, we show in Appendix C, the typical posterior probability distributions of the cosmological parameters and the parameters of EoS of the dark energy constrained by TianQin.

B. Network of TianQin and LISA
For GW observations, the power of the cosmological constraints is crucially related to how well the sky positions of the sources may be determined.With a network of multiple GW detectors working simultaneously, one can greatly improve this sky localisation error.The orbital planes of the TianQin constellation and LISA constellation are pointing in different directions; thus, joint detections can break the degeneracy between longitude and latitude of GW source, and the difference in the arrival time between the two detectors can also narrow the margin of the localisation error.However, the measurement of D L suffers from the systematics arising from weak lensing and peculiar velocities, and thus can hardly be improved even in the network observation case.
In what follows, we consider the cosmological constraints from a network of TianQin and LISA.As shown in Figs.7 and 8, compared with the TianQin alone scenario, the network of TianQin and LISA can consistently improve the constraints on the cosmological parameters.This improvement benefits from both the increased detection numbers with more detectors (as illustrated in Table I) and the better localisation capability of the network (as illustrated in Fig. 2).
Figure 7 illustrates how the relative precision of constrained cosmological parameters corresponding to different detector network configurations, MBHB population models, and weighting methods, with a similar setup as in Fig. 5. Comparing these two figures, it is easy to see that the joint detections can effectively improve the constraint precision on the Hubble constant H 0 .Similarly, the fractional dark energy density Ω Λ is still poorly constrained, except under the Q3nod model.However, different from the TianQin alone case, here we can obtain a meaningful constraints on the fractional total matter density Ω M in all cases.Again, compared with the fiducial method, the weighted method can significantly reduce the uncertainties on the Hubble constant.
When assuming a dynamical dark energy model, the joint detection of TianQin and LISA can improve the constraints on w 0 (see Fig. 8), while w a is still poorly constrained.With the weighted method scheme, for TianQin+LISA, w 0 can be constrained to a level of 19.2%, 22.3% and 7.9% in popIII, Q3d and Q3nod model, respectively; assuming a joint detection of TianQin I+II+LISA, one can better constrain the dark energy EoS parameters, w 0 can be constrained to a level of 11.3%, 18.9% and 7.5%, respectively.Especially, in the case of adopting the Q3nod model, the constraint precisions on the parameters of EoS of dark energy become comparable to those from current EM observations [24,155].Thus, the GW observations can provide an independent verification of our current understanding of the Universe [25,27,28,[156][157][158][159].
Table II summarises the relative precision on H 0 , Ω M , and Ω Λ , as well as w 0 and w a , under various detection scenarios.Each result is obtained by marginalizing over the other parameters, while all three parameters of ΛCDM model are fixed when constraining w 0 and w a .In order to mitigate the impact of random errors, we independently generate 24 GW event catalogs for each detector configuration/MBHB population model/weighting method, and present the mean values and median values of the relative constraint precisions of these five parameters.This table shows that the positional and bulge luminosity weighting of the host galaxies is very helpful for improving the precision of the parameter estimation.And a network of TianQin and LISA yield better constraints than both TianQin and LISA.In the most ideal scenario, H 0 , Ω M , and Ω Λ can be estimated with a relative precision of 1.7%, 12.4%, and 17.6%, while w 0 and w a can be constrained to 7.5% and 56.7%, respectively.

C. EM-bright scenario
Some literature has suggested that MBHB mergers can be accompanied by x-ray, optical, or radio activity [83,160], although it is not yet certain about the physical properties of the EM counterpart for GW events.In most cases, the gas-rich environment near the MBHBs is considered to be responsible for such EM transients.For example, in the late stage of binary evolution, the gas within the binary orbit would be driven inward by the inspiralling MBHB.The accretion rate can exceed the Eddington limit, forming high-velocity outflows and emitting strong EM radiation [95].If the component MBHs are highly spinning with aligned spin, the external disk can extract energy from the orbiting MBHB until merged, forming dual jets and observable emissions in a way similar to the Blandford-Znajek mechanism [96].General relativistic magnetohydrodynamic simulations of magnetized plasma show that MBHBs can amplify magnetic fields by strong accretion and emit strong EM signals [161].Some have argued that jets are increasingly dominated by magnetic fields, and can lead to transient emission after the merger [97].Hydrodynamic simulations including viscosity indicate that, before merger, the MBHB will transfer orbital energy into internal shocks and emit through x-ray radiation [98].
In view of these literatures, it is interesting to consider the optimistic scenario, where an EM counterpart of each MBHB merger is identified, and these counterparts are further used to extract redshift information from the host galaxy observation [83].This possibility would be enhanced by a number of current or planned EM facilities with large field of view and high sensitivity, including Einstein Probe [162,163], Chinese Space Station Telescope (CSST) [145], Euclid telescope [146], Vera Rubin Observatory [164], Five-hundred-meter Aperture Spherical radio Telescope (FAST) [165], and Square Kilometre Array (SKA) [166].We remark also that this bright standard siren analysis can serve as a lower limit for the precision of cosmological parameter estimation in the dark standard siren scenario -at least in the absence of any additional mitigating strategy to deal with the impact of weak lensing and/or peculiar velocity errors.In our bright standard siren analysis, we apply the same selection criteria as [83], and only consider the GW events with ρ 8 and ∆Ω 10 deg 2 (and z < 3, additionally).
In Fig. 9, we illustrate the expected probability distribution of (H 0 , Ω M , Ω Λ ), while in Fig. 10 we show the expected results for the CPL parameters (w 0 , w a ).In both cases, we consider three detector configurations, namely TianQin, TianQin I+II, and TianQin+LISA, represented by the grey shadow, black line, and magenta shades, respectively.We also consider the popIII, Q3d, and Q3nod models.Correspondingly, the marginalised one dimensional relative errors of the parameters for various detector configuration are listed in Table III.
As expected, with the identified EM counterparts, one can significantly improve the ability of the standard sirens to constrain the cosmological parameters.For TianQin, the relative precision of H 0 can be as small as 1.9%, which translates into an absolute precision of about ∆H 0 = 1.3 km/s/Mpc.For TianQin I+II, the relative precision of H 0 can be as small as 1.4%, which exceeds the constraint accuracy for TianQin+LISA in the dark standard siren scenario.If TianQin I+II+LISA is implemented, the relative precision of H 0 can be as small as 1.2%, which translates into an absolute precision of about ∆H 0 = 0.8 km/s/Mpc.Furthermore, the fractional density parameter Ω Λ and the EoS parameter w a for the dark energy, which in the dark standard siren scenario are hardly constrained, can also be constrained to a relative precision of as small as 9.5% and 33.0% when the EM counterpart available.These constraints are comparable to those from analysis of other cosmological EM survey data [155,167].For a given GW event, its ability to constrain the cosmological parameters depends on two factors: its spatial localization accuracy, and the cosmological evolution between the source and the observer.
0.1 0.3 0.5 0.7 0.9 1.1 1.3 1.5 1.7 1.9 2.1 2.3 2.5 2.7 2.9 In Fig. 11, we show the distribution of sky localization errors for MBHBs mergers at different redshifts.We see the sky localization area ∆Ω tends to increase as the redshift increases.Meanwhile, nearby events are usually accompanied by larger SNR, which leads to better determination of the luminosity distance.To sum up, for nearby events, the smaller localization area together with the better distance estimation leads to a smaller error box, and therefore a smaller number of candidate host galaxies.
On the other hand, given the same ∆D L and ∆Ω, the GW signals at different redshifts lead to different constraints on the cosmological parameters.In Fig. 12, we illustrate the degeneracy of the cosmological parameters by showing the evolution of correlation coefficients between different pairs of parameters.To do this, we estimate cosmological parameters with a sample of EM-bright GW events, assuming a relative error on luminosity distance σ D L /D L = 0.1.
We find that, the H 0 constraint from high redshifts suffers from strong degeneracy with Ω M and Ω Λ , and this degeneracy becomes less significant only at very low redshifts.Notice that the correlation coefficient between Ω M and Ω Λ switches its sign at z ∼ 3.5, suggesting that observations of both higher and lower redshift events are needed to alleviate the degeneracy.Also, the correlation coefficient between w 0 and w a is approaching total anticorrelation at z 1.These two parameters cannot be precisely measured, partly because of the strong anticorrelation between them.
In conclusion, a small number of low redshift GW sources is sufficient to provide a good constraint on the H 0 , thanks to the precise position estimation, as well as the weaker parameter degeneracies.Meanwhile, high redshift events are useful to better constrain other cosmological parameters like Ω M and Ω Λ .For cases where the EM counterpart can be uniquely identified, the GW events with large redshift can give tight constraints on the cosmological parameters; otherwise the large positional error would weaken such constraints.

B. Effects of redshift limit
Throughout this study we have applied a cutoff for events beyond redshift z = 3.To some extent, this cutoff reflects the incompleteness of galaxy catalogs at higher redshifts.It is very challenging to pursue relatively complete galaxy catalog at high redshift.Such as, sky surveys like the Sloan Digital Sky Survey (SDSS) [168][169][170] and Dark Energy Survey (DES) [171] can reliably map galaxies as far as z 1.2.
Our paper highlights the need for more detailed galaxy surveys with redshift limit of z > 3. Actually, it is not beyond the imagination that, once MBHB mergers are routinely detected, intensive observations would be triggered to map the higher redshift Universe -at least within the sky localization regions of the observed mergers -and thus provide a more complete catalog of galaxies therein [172].
Moreover, the redshift limit for quasar observations is much larger than that for galaxies.For example, the redshift limit of the quasar catalog mapped by SDSS can reach z 4 [169,170].Meanwhile, some have argued that quasars can host the MBHBs [173][174][175].In addition to the quasars, one can also explore the potential of using Lyman-α forest effect to obtain redshift information [176][177][178][179].It is still possible to obtain statistical redshift for the high-redshift GW sources.

C. Importance of bulge luminosity information
Throughout this paper, we have incorporated both position and bulge luminosity information for the calculation of weights that we assign to candidate host galaxies.Here we demonstrate the importance of the bulge luminosity alone, by weighting only on spatial location information.As shown in Fig. 13, compared with the fiducial method, the relative error on H 0 shrinks when sky location information is included, but the bulge luminosity weight consistently improves the estimation precision.For mergers with higher redshift, the importance of the bulge luminosity is weakened, due to the large localization error volume and the incompleteness of the galaxy catalog.The large number of galaxies within this volume, combined with large uncertainty associated with the bulge luminosity, decreases the effectiveness of the weighting scheme.On the other hand, the galaxy catalog is incomplete at high redshift, and the introduction of the luminosity function and the bulge luminosity function only compensates the weights of high redshift bright galaxies to some extent [see Eqs.(30) and ( 31)], which does not help to weight the "correct" host galaxies for the high redshift GW sources with lower mass.
,QGH[RIVLPXODWLRQ Comparison of the constraints on the Hubble constant H0 by an individual nearby event for the following scenarios: with EM counterpart (blue star), using the position+luminosity weighted method (red square), and using position-only weighted method (cyan triangle).Each index represents one nearby GW event.The green line represents the H0 error measured by P lanck using CMB anisotropies [24], and the gray line represents the H0 error measured by the SH0ES project using type Ia supernovae data [28].These events are selected from the sample pool of the Q3nod model with the conditions of both z 0.5 and ∆Ω 5 × 10 −3 deg 2 .
However, the bulge luminosity weight plays a crucial role especially when nearby events are considered, since only a relatively small number of galaxies is located within the error box.In this case, therefore one host galaxy can be largely identified through the M MBH − L bulge relation.Fig. 14 illustrates the constraints on the Hubble constant for the nearby events with both z 0.5 and ∆Ω 5 × 10 −3 deg 2 in three scenarios, i.e., with EM counterpart, using the position plus luminosity weighted method, and using position-only weighted method.We can observe that, by using the M MBH − L bulge relation, the constraint precision on the Hubble constant in the dark siren scenario can be greatly improved, sometimes even comparable to scenarios with EM counterpart.Since we also expect that it would be easier to obtain more complete information on the bulge luminosity for nearby galaxies, this fact highlights the importance and potential efficacy of the bulge luminosity information.

D. Scope of the MMBH − L bulge relation and the luminosity function
We have demonstrated that performing a weighting process of candidate host galaxies according to the M MBH − L bulge relation is vital to improve the precision of the measured cosmological parameters.In this paper, a default assumption is adopted in our calculation, which is that the M MBH − L bulge relation applies to all MBHBs independent of the mass of the GW source.This relation behaves differently at the low-mass end than at the high-mass end [103,104], and we have used different values of intrinsic scatter at the high and low-mass ends, see Eq. (28).In this part we discuss the scope of application of this relation.
For TianQin, the MBHB events with M > 5 × 10 6 M account for about 18%, 16% and 21% of the total detectable events under popIII, Q3d and Q3nod models, respectively; for LISA, these percentages are about 44%, 36% and 40%, respectively.The high-mass MBHB GW sources account for quite a large proportion of the total detectable events.For the low-mass MBHs with M < 10 6 M , their host galaxies show no co-evolution with the central MBHs [104,180].The low-mass MBHs are more likely to be remnants of black hole seeds, for a given MBH mass, the possible bulge luminosity spans a larger range [104].
Bulges at the centers of galaxies can be classified as either classical bulges or pseudo-bulges, depending on whether they contain disk structures (pseudo-bulges) or not (classical bulges).Based on numerical simulations of galaxy collisions, it is generally accepted that classical bulges are made in major mergers of galaxies, whereas pseudobulges are the result of internal evolution of galaxy disks since their formation [103].It is noteworthy that the host galaxies of the low-mass MBHs that deviate from the M MBH − L bulge relation almost all have pseudobulges in their centers [104].This implies that the deviations from the M MBH − L bulge relation at the low-mass end do not unduly affect the scope of the weighted method proposed in this paper, and it is sufficient to account for this effect by adopting a larger intrinsic scatter, see Eq. (28).
In addition to the bulge luminosity information of galaxies, other observable information about galaxies such as stellar velocity dispersion and galaxy morphology may provide additional constraints on MBH mass for the lighter ones.The relation between MBH mass and stellar velocity dispersion (the M MBH − σ relation) at the low-mass end is consistent with that at the high-mass end [181].Even such information is not available a prior, telescopes would be motivated to perform deep observations after the GW detection, especially when the sky location can be very precisely determined.
We also examine other aspects of this relation.Throughout this paper we have assumed a universal M MBH − L bulge relation, obtained from fitting the observed data at z ≈ 0 [102,103].However, due to the lack of reliable models, we ignore the possible redshift evolution of this M MBH − L bulge relation, which might cause some bias in our analysis.In the weighting process, we choose a moderately large σ log 10 M so that such bias can be partially absorbed.
Another possible source of bias is introduced by ignoring the redshift evolution of the luminosity function.We manually augmented our galaxy samples, according to the modelled luminosity function, to remove the Malmquist bias.However, the luminosity function is only complete for nearby redshifts, and becomes more and more incomplete as the redshift increases [113,148,182].We therefore highlight also the need for a more thorough and precise luminosity function evolution model, which can be obtained by existing [114,183] and future ultradeep-field observations [145,146].
Notice that most galaxy catalogs constructed from cosmological surveys do not list the galactic bulge luminosity of each source.In order to use bulge luminosity to improve our weighting scheme, we need to extract this information from data.Since the luminosity of a galaxy decreases with its radius, which can be described empirically using the Sérsic function (also known as the r 1/n law) [184][185][186][187][188], we can approximately estimate the bulge luminosity from the total luminosity and the morphology of a galaxy based on the Sérsic function, but this conversion would bring additional large uncertainties.

E. Consistency check
The possibility exists that there will be only a small number of GW events available throughout the observation time, and fluctuations resulting from small number statistics could potentially bias the estimated cosmological parameters.This bias can happen especially when the true host galaxy is relatively dim, so that the galaxy survey could mistakenly associate the GW event with some other galaxy.Following [87], we discuss the impact of this bias, as well as how to perform a consistency check, which can identify and remove it.
Without loss of generality, in the following analysis we focus on the results for H 0 and w 0 .For a given simulated catalog, we first perform a cosmological parameter estimation using all GW events.Then, event by event, we remove one event (hereafter the kth event), and perform the cosmological analysis using the events, which remain, with p( Ω|I) i =k p(D i , S| Ω, I).We find that, the mis-association of the host galaxy could occur, and could severely bias our estimation.
An example of this bias is illustrated in Fig. 15, where the mis-association of one event occurs, so that posteriors on H 0 and w 0 obtained from all subcatalogs that contain this event (indicated by the grey lines) yield estimates of the parameters that deviate significantly from their true values.Notice, however, that the posteriors obtained from the subset that excludes this particular mis-associated event (indicated by the red line) show in each case an obvious deviation from the remaining posteriors, and thus provide a clear diagnostic for the potential mis-association bias.
The above analysis therefore provides a consistency check for the mis-association.If the number of events is too large to carry out such a consistency check, one can remove more events in every iteration in order to enhance its efficiency.(Moreover, we also note that when the total number of events is larger, then the impact of mis-association of a single event will also be less severe.) In each panel, the black line represents the result for the whole catalog, the grey lines represent the results for subsets that include the mis-association, and the red line represents the result for the subset that exclude the mis-association.The obvious deviation of posterior modes between the red line and the grey lines indicates the possible occurrence of a mis-association bias.

VI. CONCLUSIONS AND OUTLOOK
In this paper, we develop a Bayesian analysis framework for constraining the cosmological parameters, using simulated observations of the MBHBs mergers from space-borne GW observatories like TianQin and LISA.We obtain the luminosity distance information directly from the GW observations, while a statistical analysis of simulated galaxy catalogs from EM surveys is used to obtain the corresponding redshift information.With the identification of an explicit EM counterpart, one can indeed perform very precise cosmological measurement.However, we also show that one can still obtain useful cosmological constraints even in the dark standard siren scenario where no EM counterpart exists -provided that Malmquist bias is properly accounted for.Furthermore, if we include the localization and mass information from the GW observation to inform the weighting of candidate host galaxies, making use of the relation between the central MBH mass and the bulge luminosity of the host galaxy, we can significantly improve the cosmological constraints.
For the dark standard siren scenario, we consider two weighting schemes, namely the fiducial method, one with uniform weights for all galaxies within the error box, and the weighted method, the weights related to its location and bulge luminosity.With the weighted method scheme, the precision of cosmological parameter estimates can be greatly improved.In this scheme, TianQin can constrain the Hubble constant H 0 to a precision of 6.9%, 6.5%, and 3.3%, and the CPL parameter w 0 to a precision of 36.7%, 27.2%, and 13.8%, for the popIII, Q3d and Q3nod MBHB population models respectively.For TianQin I+II, the H 0 precision can be improved to 6.0%, 6.0%, and 2.0%, and w 0 precision can reach 26.6%, 29.6%, and 8.1% -again for popIII, Q3d, and Q3nod models, respectively.The other parameters like Ω M and Ω Λ need a larger number of events to be significantly constrained.Under the Q3nod model, using the weighted method, Ω M and Ω Λ can be constrained to an accuracy of 34.9% and 26.1% for TianQin, and 25.1% and 21.5% for TianQin I+II, respectively.However, w a is always difficult to be constrained in all cases.
LISA can perform similarly to TianQin I+II, but the joint detection of TianQin and LISA can significantly improve the precision of the cosmological parameter estimates relative to the results obtained from an individual detector.Using the weighted method, for TianQin+LISA, the precision of H 0 improves to 4.7%, 5.2%, and 1.8%, and the precision of w 0 improves to 19.2%, 22.3%, and 7.9% for popIII, Q3d, and Q3nod models, respectively; for TianQin I+II+LISA, the precision of H 0 improves to 3.4%, 4.7%, and 1.7%, and the precision of w 0 improves to 11.3%, 18.9%, and 7.5%, respectively.
For the EM-bright standard siren scenario, the identification of the EM counterpart can help to pinpoint the redshift, and one can not only gain tighter constraints on H 0 and w 0 , but also obtain meaningful constraints on the other cosmological parameters, including Ω M , Ω Λ , and w a .
The constraints on the Hubble constant H 0 are mainly derived from a few events at low redshift, but high redshift GW events play an important role of smoothing out fluctuations in the H 0 posterior, as well as helping to constrain the values of Ω M and Ω Λ .The CPL model describes the evolution of the dark energy EoS with redshift, so a combination of the GW events at different redshift is needed to constrain w 0 and w a .We also discuss the application of consistency checks on the estimated cosmological parameters, so that the potential bias caused by a single mis-association of a GW event with its host galaxy could be identified.
There are a number of assumptions, which we adopt that could affect the applicability of our method.For example, we depend strongly on the availability of a relatively complete galaxy catalog, reaching out to a high redshift.While in reality, obtaining a relatively complete galaxy catalog at z > 1.5 is very challenging, when one considers the small localization area that will be provided by TianQin and/or LISA it seems likely that deep-drilled, a more complete galaxy surveys triggered by GW detections could be carried out in order to alleviate incompleteness issues.We also assume a simple relationship between the MBHs mass and the bulge luminosity.Knowledge of the redshift evolution of this M MBH −L bulge relation could also improve the cosmological constraints by reducing any potential bias arising from adopting this simple relationship.Extra information like stellar velocity dispersion can be used to improve the estimation of MBH mass at the low-mass end through the M MBH − σ relation.
In the future there is scope to extend our paper in several ways.For example, a new population model of MBHB mergers can be added for the analysis [189].The possible existence of strong gravitational lensing events can also help with pinpointing the redshift of GW events [190,191].We also plan to extend the GW cosmology study to more types of GW sources [67,68].and Bai-Tian Tang for helpful discussions.

FIG. 2 .
FIG.2.Spatial localisation error distribution for TianQin (solid line) and the LISA-TianQin network (dot-dashed line).The top panels show the error due to GW observations, with the panel (a) showing the sky localisation error, and the panel (b) showing the relative uncertainty on luminosity distance.In the bottom panels, we illustrate uncertainties of luminosity distances, also considering uncertainties due to weak lensing and peculiar velocity.The panel (c) shows the uncertainty evolution over redshift summarised from 1000 simulations, and the panel (d) shows similar information as the panel (b) but with all uncertainty sources considered.The results with TianQin I+II are quite similar to those of TianQin so for clarity we do not include them here.

FIG. 3 .
FIG.3.An example spatial localization error for a GW source.The cyan volume represents the spatial error due to the GW observation, while for the orange volume, we convert the [zmin, zmax] back into luminosity distance assuming the correct cosmology.The gray shadows are the projections of the spatial localization error, the red dot is the observer, and the purple star is the true input position of the GW source.

FIG. 4 .
FIG.4.Example: Redshift distribution of galaxies for an example GW source.Top panel: Considering no debias correction, and assuming all galaxies within the error box as equally likely to be the host galaxy.Middle panel: Same as the top panel but including a correction for Malmquist bias.Bottom panel: Corrected for Malmquist bias, and assigning galaxy weights according to their position as well as luminosity.The vertical orange dashed line represents the real redshift of the GW event, and the horizontal blue dashed line represents the prior distribution where galaxies have constant density over comoving volume.

FIG. 9 .
FIG.9.Typical results for the estimation of cosmological parameters under the optimistic scenario with an explicit EM counterpart.Contours show the 68.27% and 95.45% confidence levels, assuming different detector configurations of TianQin (grey shadow), TianQin I+II (black line), and TianQin+LISA (magenta shade), respectively.The three subplots correspond to adopting popIII, Q3d, and Q3nod respectively as the underlying model for MBHB mergers.

FIG. 10 .
FIG.10.Typical results for the estimation of dark energy EoS parameters under the optimistic scenario with an explicit EM counterpart.Contours show the 68.27% and 95.45% confidence levels, assuming different detector configurations of TianQin (grey shadow), TianQin I+II (black line), and TianQin+LISA (magenta shade), respectively.The three subplots correspond to adopting popIII, Q3d, and Q3nod respectively as the underlying model for MBHB mergers.The other cosmological parameters, H0, ΩM , and ΩΛ are fixed.

FiducialFIG. 13 .
FIG. 13.Comparison of the relative errors on the Hubble constant H0 under different weighting methods, assuming GW detections with TianQin I+II under the Q3nod MBHB model and no EM counterpart.Left panel: Each point represents the relative error on H0 for an individual simulation corresponding to the fiducial (black dot), the position-only weighted (cyan triangle), and the position+luminosity weighted (red square) method, respectively.Right panel: Histogram of the relative error on H0 for the three different weighting methods.

FIG. 15 .
FIG. 15.A consistency check on an example catalog, showing marginalized posterior distributions of H0 [panel (a)] and w0 [panel (b)].In each panel, the black line represents the result for the whole catalog, the grey lines represent the results for subsets that include the mis-association, and the red line represents the result for the subset that exclude the mis-association.The obvious deviation of posterior modes between the red line and the grey lines indicates the possible occurrence of a mis-association bias.
FIG. C16.Typical corner plots of the posteriors for the parameters h, ΩM , and ΩΛ constrained by the detections of TianQin, comparing the fiducial method (left column) and the weighted method (right column).The top, middle, and bottom rows correspond to adopting popIII, Q3d, and Q3nod respectively as the underlying model for MBHB mergers.In each subplot, the three panels at the lower left show the two-dimensional joint posterior probabilities of h − ΩM , h − ΩΛ, and ΩM −ΩΛ, with the contours represent confidence levels of 1σ(68.27%)and 2σ(95.45%),respectively; the upper, middle and right panels show the one-dimensional posterior probabilities of the corresponding parameters, after marginalization over the other parameters, with the dashed lines indicate 1σ credible interval.In each panel the solid cyan lines mark the true values of the parameters.

TABLE I .
Expected total detection rate of GW events with z < 3 and SNR ρ > 8, detected over the entire observation time based on the MBHB population models for five different detector configurations: TianQin, TianQin I+II, LISA, TianQin+LISA, and TianQin I+II+LISA.

TABLE II .
Expected relative precision on (H0, ΩM , ΩΛ) and (w0, wa) constraints for the three MBHB models, assuming different configurations of detectors or networks.When constraining H0, ΩM , and ΩΛ, no evolution of dark energy is assumed, while all other cosmological parameters are fixed when studying the CPL parameters.Numbers are shown for mean values (median values in brackets).The symbol "−" indicates that the corresponding parameter is not effectively constrained under the corresponding condition.

TABLE III .
Constraints on the cosmological parameters and dark energy EoS parameters under the optimistic scenario that an EM counterpart is observed, considering different detector configurations and underlying MBHB merger models.Notice that the cosmological parameters (H0, ΩM , ΩΛ) and CPL parameters (w0, wa) are studied separately.Each result is obtained based on 108 replicated and independent random realizations, the error represents 68.3% confidence interval.