Measuring the Hubble Constant of Binary Neutron Star and Neutron Star–Black Hole Coalescences: Bright Sirens and Dark Sirens

Observations of gravitational waves (GW) provide us with a new probe to study the Universe. GW events can be used as standard sirens if their redshifts are measured. Normally, standard sirens can be divided into bright/dark sirens according to whether the redshifts are measured by electromagnetic (EM) counterpart observations. First, we investigate the capability of the 2.5 m Wide-Field Survey Telescope (WFST) to take follow-up observations of kilonova counterparts. For binary neutron star (BNS) bright sirens, WFST is expected to observe 10–20 kilonovae per year in the second-generation GW detection era. As for neutron star–black hole (NSBH) mergers, when a BH spin is extremely high and the neutron star (NS) is stiff, the observation rate is ∼10 per year. Combining optical and GW observations, the bright sirens are expected to constrain the Hubble constant H 0 to ∼2.8% in five years of observations. As for dark sirens, the tidal effects of NSs during merging provide us with a cosmological model-independent approach to measure the redshifts of GW sources. Then we investigate the applications of tidal effects in redshift measurements. We find in the third generation era, the host galaxy groups of around 45% BNS mergers at z < 0.1 can be identified through this method, if the equation of state is ms1, which is roughly equivalent to the results from luminosity distant constraints. Therefore, tidal effect observations provide a reliable and cosmological model-independent method of identifying BNS mergers’ host galaxy groups. Using this method, the BNS/NSBH dark sirens can constrain H 0 to 0.2%/0.3% over a five-year observation period.

In May 2023, LIGO/VIRGO and the Kamioka Gravitational Wave Detector (KAGRA) will start the fourth observing (O4) run with a higher detection sensitivity (Abbott et al. 2018b).And in 2025, the fifth observing (O5) will begin, with the A+ upgrade and the addition of an intrument in In-dia.At the same time, there will be a variety of advanced telescopes, such as the Large Synoptic Survey Telescope (LSST; Ivezić et al. 2019), in operation.These will significantly improve the detection capability of kilonova and broaden the application of bright sirens in cosmology.Among them, the 2.5-meter Wide-Field Survey Telescope (WFST), which is located at the Saishiteng Mountain in the Lenghu area (93 • 53 ′ E, 38 • 36 ′ N), is scheduled to be installed in the early of 2023.It has a field-of-view (FOV) of 6.55 deg 2 , 5σ limiting magnitudes of (22.40, 23.35, 22.95, 22.59, 21.64, 22.96) in the (u, g, r, i, z, w) for a 30s exposure time (Liu et al. 2023), and a median seeing of 0.75 arcseconds in the Lenghu site (Deng et al. 2021).With these parameters, WFST is expected to be one of the major follow-up telescopes in the northern hemisphere during the O4 and O5 runs of the advanced LIGO, advanced Virgo and KAGRA.In the first half of this paper, we will discuss the capability of WFST to perform follow-up observations during the O4/O5 runs and foresee the application of the observed BNS/NSBH bright sirens in the Hubble constant constraints.
Multi-messenger observations are a very direct and effective method of measuring redshift.However, among the other GW events observed by LIGO/VIRGO so far, only GW190521/ZTF19abanrhr (Abbott et al. 2020c;Graham et al. 2020) is a probable multi-messenger event.For those 'dark standard sirens ', MacLeod & Hogan (2008) and Del Pozzo (2012) proposed that their redshifts could be statistically obtained from the comparison of GW events' localization areas and the survey data.The ability of localizing GW events has been improved significantly with the operation of VIRGO, and this method has been widely discussed in recent years (Nishizawa 2017;Nair et al. 2018;Chen et al. 2018;Fishbach et al. 2019;Soares-Santos et al. 2019;Palmese et al. 2019;Yu et al. 2020;Gray et al. 2020;Palmese et al. 2020;Vasylyev & Filippenko 2020;Borhanian et al. 2020;Zhao et al. 2020;Jin et al. 2020Jin et al. , 2022aJin et al. ,b, 2023;;Abbott et al. 2021b;Finke et al. 2021;The LIGO Scientific Collaboration et al. 2021;Palmese et al. 2021;Wang et al. 2022;Song et al. 2022).In particular, this method will be most effective when the number of probable host galaxies, groups or clusters in the localization region is the only one.However, since the localization regions from GW detectors are related to d L , in order to compare them with the survey catalog, these localization regions need to be converted to redshift coordinates under priors of specific cosmological models.Therefore, due to the indeterminacy in the d L − z relation, the existence of other galaxies/groups/clusters with different redshifts in the same sky region will have a significant impact on redshift statistics.
In order to overcome the reliance on cosmological models, a natural way is to measure redshifts directly from GW data.In GW waveforms, the redshift of a binary system is usually tightly coupled to its physical mass m phy .If m phy could be measured, the coupling would break down and the redshift could also be measured.Neutron stars (NSs) under the effect of strong tidal fields will have a tidal deformation, which will have an effect on the waveform (Vines et al. 2011;Bini et al. 2012;Vines & Flanagan 2013;Abdelsalhin et al. 2018;Landry 2018;Banihashemi & Vines 2020).Since the tidal deformability λ is related to the equation of state (EOS) as well as the physical mass of the NS, if the EOS is known, then the physical mass of the NS can be measured by tidal effect and thus the redshift can be determined (Messenger & Read 2012;Del Pozzo et al. 2017;Wang et al. 2020;Chatterjee et al. 2021).
Currently, GW170817 has given a loose constraint on λ (Abbott et al. 2017c) and excluded some of the EOSs (Abbott et al. 2018a).In the future, with the third generation (3G) of GW detectors CE (Dwyer et al. 2015) and ET (Punturo et al. 2010) in operation, the accuracy of λ measurements will improve significantly.Meanwhile, observations of EM counterparts such as kilonova (Nicholl et al. 2021) and X-ray plateau (Lasky et al. 2014;Gao et al. 2016;Li et al. 2016;Ren et al. 2020;Lucca & Sagunski 2020;Bauswein et al. 2020;Lü et al. 2021), and the search for massive NSs in the Milky Way (Antoniadis et al. 2013;Alsing et al. 2018;Cromartie et al. 2020), will impose strong constraints on the NS EOS.Thus, in the 3G era, tidal effect observations will be a powerful and cosmologically model-independent method of measuring the dark siren redshifts (Wang et al. 2020).
In our previous work (Yu et al. 2020), we have discussed the ability of the GW detector network to identify stellar mass binary black hole (SBBH) mergers' host galaxy groups by comparing d L localization results with the SDSS DR7 halobased group catalog (Yang et al. 2005(Yang et al. , 2007;;Abazajian et al. 2009), and estimated the potential of these SBBH mergers to constrain the Hubble constant.Compared to galaxies, groups represent a larger physical structure and are therefore more easily identified.In addition, the use of group catalog can partially overcome the influence of peculiar velocity of the galaxies.For these two redshift measurement methods, we follow the calculations in Yu et al. (2020) and compare their ability of identifying the BNS/neutron star-black hole (NSBH) mergers' host groups and constraining the Hubble constant.To exclude the influence of external factors such as edge effect and incompleteness effect, we imply a mock galaxy group catalog based on the Uchuu simulation (Ishiyama et al. 2021) instead in this work, which covers the whole sky and has a much higher upper redshift limit.
Our paper is organized as follows.In Section 2, we introduce the GW waveform, the detector response and the Fisher information matrix method.We introduce the kilonova model used in the analysis of bright sirens in Section 3 and then estimate the application of bright sirens in the Hubble constant constraints in Section 4. In Sections 5 and 6, we analyze the dark sirens and their cosmological applications.And in Section 7, we summarize our results.

The detector's response
In the paper, we mainly consider the array with N d GW detectors, whose spatial locations are denoted as r I with I = 1, 2, • • • , N d .For an incoming GW signal propagating along n, the response of the i-th detector in the frequency domain can be written as where h + (t) and h × (t) are two wave polarizations in the transverse traceless gauge, F + I and F × I are two beam-pattern functions, which depend on the source's right ascension α, declination, δ, the polarization angle, ψ, the position and the orientation of the interferometer's arms.In this paper, for bright sirens, we discuss the network of LIGO (Livingston), LIGO (Handford) and Virgo, and the network of five 2G GW detectors with the addition of KAGRA and LIGO-India (Unnikrishnan 2013).We denote these two networks as LHV and LHVKI, respectively.For dark sirens, we focus on a network of three 3G GW detectors, which consists of ET, CE and an assumed CE-type detector located in Australia.Their parameters and noise curves can be found in our pervious paper, Yu et al. (2020) 1 .For ET, we use the proposed ET-D project noise curves (Punturo et al. 2010).For CE and assumed CE-type detector, we consider the proposed noise curve in Abbott et al. (2017f) and Dwyer et al. (2015).We denote this network as CE2ET.
In this paper, we adopt the restricted post-Newtonian (PN) waveform for the nonspinning systems in an inspiralling stage (Cutler et al. 1993;Damour et al. 2001;Itoh et al. 2001;Blanchet 2002;Blanchet et al. 2002Blanchet et al. , 2004a,b;,b;Itoh & Futamase 2003;Itoh 2004a,b;Blanchet & Iyer 2005;Sathyaprakash & Schutz 2009).For a binary system with component masses m 1 and m 2 , total mass M = m 1 + m 2 , mass ratio η ≡ m 1 m 2 /M 2 , chirp mass M c ≡ Mη 3/5 , inclination angle ι and luminosity distance d L , the waveform is given by where t c is the merging time.The amplitude A I is given by where M c ≡ (1 + z)M c is the "redshift chirp mass".Meanwhile, after defining M = (1 + z)M, the functions φ and ϕ I,(2,0) are given by where we use the 3.5 PN approximation for the phase, and φ c is the merging phase.The parameters φ i can be found in Sathyaprakash & Schutz (2009).For a BNS merger, the tidal deformation term Φ tidal ( f ) is given by (Vines et al. 2011; Messenger & Read 2012) where χ a ≡ m a /(m 1 + m 2 ), and λ a (m a ) characterizes the changes of the induced quadrupole with an external tidal field.In this paper, we use the λ − m relation under the linear approximation mentioned in Wang et al. (2020), where B and C are two parameters determined by the EOS of NSs.Following Wang et al. (2020), we discuss a sam-ple of 13 NS EOSs,alf2,ap3,ap4,bsk21,eng,H4,mpa1,ms1,ms1b,qmf40,qmf60,sly,wff2,and 4 quark star EOSs, MIT2, MIT2cfl, pQCD800, sqm3 in this paper ( Özel & Freire 2016;Zhu et al. 2018;Zhou et al. 2018;Xia et al. 2021).The corresponding fitted values of B and C can be found in Table 2 of Wang et al. (2020).The maximum stable NS mass, M TOV (Tolman-Oppenheimer-Volkoff mass) and the radius of a 1.4 M ⊙ NS, R 1.4 can also be found in this paper.For a NSBH merger, the tidal deformation term Φ tidal ( f ) can be obtained by simply setting λ 1 = 0 in Eq. ( 6), where a = 1 represents the BH component.
From the PN waveforms above, we can see that in all terms except Φ tidal ( f ) , redshift z and mass M are tightly coupled through the redshift mass M = (1 + z)M.With the addition of the tidal deformation term, the coupling between redshift and mass will break down, and we can obtain the redshift of GW source from the phase observations.Therefore, for a BNS/NSBH merger, the response of the interferometer depends on twelve parameters, θ = {α, δ, ψ, ι, M, η, t c , φ c , log(d L ), z, B, C}.Using the method in Wen & Chen (2010), we can obtain the Fisher matrix Γ i j and covariance matrix (Γ −1 ) i j for these 12 parameters.The signalto-noise ratio (SNR) of the GW signal can also be obtained from the GW detector's response (see Zhao & Wen 2018 for details of the Fisher matrix and SNR).In this work, for dark sirens we choose SNR > 12 as the threshold for GW detection.For the bright sirens, we relax this threshold to SNR > 8 due to the existence of EM counterparts.

Samples of GW events
In this paper, we adopt the BNS mergers' redshift distribution model in our previous work Yu et al. (2021), with the star formation rate (SFR) model from Yüksel et al. (2008), the log-normal time delay model from Wanderman & Piran (2015), and the local merger rates updated to GWTC-3, R BNSmergers,0 = 13 − 1900 Gpc −3 yr −1 , and R NSBHmergers,0 = 7.4 − 320 Gpc −3 yr −1 (The LIGO Scientific Collaboration et al. 2021b).Since there is a large difference between the upper and lower bounds of the merger rates of GWTC-3, we use the geometric mean of the upper and lower bounds, R = R low R high in our calculations to simplify the analysis, which are R BNSmergers,0 = 157.2Gpc −3 yr −1 and R NSBHmergers,0 = 48.7 Gpc −3 yr −1 .After adopting these two local merger rates, the total numbers of z < 0.2 BNS and NSBH mergers per year in the whole sky N yr are about 600 and 180, respectively.These two numbers are roughly in line with the predictions in Iacovelli et al. (2022).
We assume a uniform distribution of NS masses between 1.2 M ⊙ and 2.0 M ⊙ .This assumption is consistent with the observational constraints of GWTC-3 (The LIGO Scientific Collaboration et al. 2021b), and keeps the mass of the sample below the maximum stable NS mass M TOV of the EOS mentioned in Wang et al. (2020) (except for sqm3, M TOV = 1.99 M ⊙ ).For BH, we follow The LIGO Scientific Collaboration et al. (2021b) and use a power-law mass distribution model in the interval [2.5 M ⊙ , 50 M ⊙ ] with an index of -3.4,supplemented by a Gaussian peak at 34 +4.0 −3.0 M ⊙ .The distribution of inclination angles ι is ∝ sin ι.All of the BNS and NSBH mergers in our simulation are nonspinning If no special instructions.

KILONOVA MODEL
During the merger of BNS and part of NSBH binaries, neutron rich ejecta will be ejected through tidal interactions (Rosswog et al. 1999;Sekiguchi et al. 2015), material squeezed at the contact interface (Oechslin et al. 2007;Bauswein et al. 2013;Hotokezaka et al. 2013) and disk outflows (Surman et al. 2008;Metzger et al. 2008;Nedora et al. 2019).Rapid neutron capture (r-process) nucleosynthesis will happen in these ejecta and produce heavy elements (Symbalisty & Schramm 1982) and the decay of these heavy elements will power a rapidly evolving, roughly isotropic thermal transient named as 'kilonova' (Metzger et al. 2010).Since kilonovae have a much larger viewing angle compared to GRBs, they are believed to be the main EM counterpart search targets for low redshift BNS and NSBH mergers.
POSSIS (Bulla 2019) is a three-dimensional Monte Carlo code for predicting spectra and light curves of supernovae and kilonovae.With this code, Dietrich et al. (2020) and Anand et al. (2021) provide kilonova simulation grids based on a threecomponent BNS model and a two-component NSBH model, repectively.In this section, we use POSSIS code to obtain the kilonova light curves.

Kilonova emission from BNS mergers
For the three-component BNS model in Dietrich et al. (2020), the first and second components are ejected dynamically through tidal interactions (Rosswog et al. 1999;Sekiguchi et al. 2015) and material squeezed at the contact interface (Oechslin et al. 2007;Bauswein et al. 2013;Hotokezaka et al. 2013), respectively.The former component is neutron rich and concentrated near the equatorial plane, with a high opacity and named as 'red' component.The second component has a much lower neutron abundance and opacity due to the e ± captures and neutrino irradiation (Sekiguchi et al. 2015;Radice et al. 2016), and named as 'blue' component.The third component comes from the post-merger disk outflows (Surman et al. 2008;Metzger et al. 2008;Nedora et al. 2019).Its velocity is much slower than dynamical ejecta, and its opacity lies between the first two components.After giving different values of opacities to the three components and characterizing the velocity with ejecta mass, Dietrich et al. (2020) simulate Spectral Energy Distributions (SEDs) and light curves with four free parameters, which are the mass of the dynamical ejecta M dyn , the mass of post-merger disk-wind ejecta M wind , the half-opening angle of the tidal dynamical ejecta Φ, and the viewing angle θ obs .
For M dyn , we use the fitting formula from Coughlin et al. (2019a), is the C i is the compactness of NS with radius r i .The term [1 ↔ 2] represents exchanging the subscripts.When binary masses are unequal, the less massive NS will be disrupted by the tidal forces before the collision, which will suppress the production of shocks (Bauswein et al. 2013;Hotokezaka et al. 2013;Lehner et al. 2016;Dietrich & Ujevic 2017).Here we use the simulations by Sekiguchi et al. (2016) to estimate the fractions of the red and blue components.Nicholl et al. (2021) presented a polynomial fit to these simulation results.For BNS mergers with q ≳ 0.8, the fit can be written as with fitting parameters a = 14.8609, b = −28.6148,c = 13.9597, and for mergers with q ≲ 0.8, f red ≈ 1.In the POS-SIS model, the half-opening angle of the red component ±Φ is left as a free parameter.Here we assume that the volume ratio of the red and blue components is proportional to their mass ratio.Therefore, after using the estimations of AT2017gfo, Φ = 49.50 • (Dietrich et al. 2020) as a calibration, Φ can be expressed in terms of f red .
For disk mass, we adopt the fitting formula from Dietrich et al. (2020), where M thr is the threshold mass of a NS prompt collapsing into a BH, which can be parameterised as a and b in Eq. ( 11) are given by where the best fitting values of parameters in Eq. ( 11) can be found in Dietrich et al. (2020), which are a o = −1.581,δa = −2.439,b o = −0.538,δb = −0.406,c = 0.953, d = 0.0417, β = 3.910, q trans = 0.900.Approximately 10-50% of the disc mass will be ejected by viscously-driven winds (Siegel & Metzger 2018;Radice et al. 2018), we choose a typical fraction, ϵ disk ≈ 0.2 in this work.For a BNS merger and a given NS EOS, we can then use the above equations to calculate its M dyn , M wind and Φ.With these parameters and the values of θ obs and d L , we obtain the light curves by interpolating the grids in Dietrich et al. (2020).

Kilonova emission from NSBH mergers
For NSBH systems, their BH components differ greatly from the properties of NSs during merger.So they will produce kilonova emissions that are quite different compared to BNS mergers.Firstly, NSs only produce ejecta when they are tidally disrupted outside the radius of the innermost stable circular orbit (ISCO) R ISCO , or they will be swallowed whole.The normalized ISCO radius, RISCO = R ISCO /M BH is given by (Bardeen et al. 1972 where M BH and χ BH are the mass and the dimensionless spin parameter of a BH, respectively.The tidal disrupted radius of the NS in Newtonian theory is approximately equal to (Zhu et al. 2020) where R NS and M NS are the radius and mass of a NS.It can be seen that the significant kilonova emission occurs only when the BH has a small mass and a large spin.Secondly, the dynamical ejecta are primarily produced by tidal disruption and concentrated around the equatorial plane for NSBH mergers (Kawaguchi et al. 2015), so Anand et al. (2021) consider only two components in their NSBH's kilonova model and the value of Φ is fixed at 30 • .
For the masses of dynamical ejecta from NSBH mergers, we adopt the Zhu et al. (2020)'s formulation based on 66 numerical relativistic simulation results, where C NS is the compactness of the NS, and we apply the fitting results of Gao et al. (2020) for the NS's baryonic mass M b NS in this work, , 0 The disc mass M disk can therefore be obtained by Then, similar to the BNS mergers, the kilonova light curves generated by NSBH merging can be obtained by interpolating the grids from Anand et al. (2021).

BRIGHT SIRENS
In our previous work Yu et al. (2021), we have discussed the applications of GW-GRB multi-messenger observations in constraining the EOS of dark energy.In this section, we continue the discussion of the bright sirens in cosmology.For kilonovae, which are much less luminous than GRBs but have much larger observable angles, they will be the main EM counterpart search targets for low redshift BNS and NSBH mergers.Due to limitations in the detection capabilities of GW detectors and telescopes, in this section we simulate 10000 low redshift (z < 0.2) BNS and NSBH mergers, respectively, and calculate the strength of their GW and kilonova emissions, as well as the probability of being detected by multi-messenger.Subsequently, by multiplying the normalization factor N yr /10000, we calculate the multi-messenger detection rates of the BNS and NSBH mergers and their potentials in the Hubble constant constraints.
4.1.Kilonova searching For each BNS and NSBH merger, we calculate the SNR of its GW signal, derive its covariance matrix Cov[α, δ, log(d L )] through the Fisher matrix method, and assume that its actual location is normally distributed.For each GW event with SNR > 8, we pixelate its localization area with the HEALPix pixelization algorithm (Górski et al. 2005;Zonca et al. 2019).Through the HEALPix code, the whole sky is divided into 12×nside 2 pixels, we choose nside = 512 and obtain the probability of the source lying at each pixel from the covariance matrix.
In the follow-up observations of GW events during the O2 and O3 of the advanced LIGO and advanced Virgo, many telescopes, including ZTF, DECam, etc., have designed their observation plans based on the early warning localization results of GW events, and prioritized the observation of high probability sky region to improve the efficiency of follow-up observations.In this work, we consider the effect of the telescope's observation plan, and generate a WFST observation plan for each simulated merger by combining the results of GW source localization and other properties of the event.Then we estimate the multi-messenger observation rate by combining the observation plan of each event.
To generate the observation plan, we use the Python-based package gwemopt, which was first developed by Coughlin et al. (2018) to optimize the search for EM counterparts of the GW events, and has been refined in the subsequent observations (Almualla et al. 2021).In O3, many follow-up observation plans are developed with the help of gwemopt for both single-telescope observations and joint observations with multiple telescopes (e.g., Coughlin et al. 2019b;Kasliwal et al. 2020;Antier et al. 2020;Anand et al. 2021;Frostig et al. 2022).Given a position of a telescope, gwemopt takes into account the influence from the Sun, the Moon and the Milky Way disk.In our analysis, we adopt the following settings: i) it is night when the sun's altitude is below −15 • ; ii) the maximum airmass for telescopic observations is 2.0; iii) the telescope only observes sky areas with galactic latitude |b| > 15 • .In addition, to reduce the influence of the Moon, we subtract the observation area near the Moon, the size of which is determined by the lunar phase.
gwemopt provides various algorithms in the following computational processes, i) skymap tiling, ii) time allocations and iii) scheduling.Coughlin et al. (2018) discussed in detail the efficiency of various combinations of these algorithms, among which the combination of MOC (multiorder coverage) algorithm, power law algorithm and greedy algorithm is the most efficient, so we choose this combination of algorithms in the simulations.The details of these algorithms can be found in Coughlin et al. (2018).
For each GW event, we set the observation window to three days after merging, during which the telescope will repeatedly cover the localization area.At the end of the gwemopt running time, the code will output a list with the coordinates of each exposure, the observation time, and the probability of observing a kilonova in the i-th tile P i kn .For the j-th observation in the i-th tile, its detection depth d i, j L, obs can be obtained from the absolute magnitude of the kilonova and the limiting magnitude of the telescope.Then P i kn can be obtained from where P i sky is the probability that the GW source is located at the i-th tile, CDF i (d i d L ) is the CDF of kilonova's luminosity distance at the i-th tile which can be obtained from the Fisher matrix, d i L, obs is the maximum value of d i, j L, obs .Finally, the probability of this kilonova's observation is 4.2.BNS mergers For the LHV and LHVKI networks, there are about 7.4% and 15.1% z < 0.2 BNS can be detected, respectively.In the left panels of Figure 1, we show the distributions of these observable BNS mergers' redshifts, ∆Ω and ∆ ln(d L ).For LHV, the observed BNS mergers are mainly distributed around z ∼ 0.08, and for LHVKI, z ∼ 0.11.The typical localization areas given by LHV are about 40 deg 2 (90% CL).With a longer baseline and higher detection SNR, the typical localization areas of the LHVKI array has improved to ∼10 deg 2 .Abbott et al. (2018b) predicted that during the O4, the median 90% credible region for localization area of BNS is about 33 deg 2 .The LHV results of ours are very close to them.After normalization, the detection numbers of LHV and LHVKI are 44 and 91 in one year of full-operation period, respectively.
However, a recent work, Petrov et al. (2022), estimated this value with a data-driven method.They found that with the improvement of data analysis methods in O3, many signals that were only triggered by a single detector or had a low SNR were also identified, which on the one hand increased the number of detected events, but on the other hand also decreased the average localization accuracy.In their estimation, the median 90% credible region for localization area of BNS is 1820 deg 2 for O4, which is far greater than our results as well as those of Abbott et al. (2018b).Thus our simulation may lose these 'bad-localized' events.However, on the orther hand, they also found these method improvements would basically not change the rate of the well-localized events.Considering that the EM counterparts of those bad located events are difficult to be searched, we neglect their effects in this work and leave them for future work.
For each BNS merger with SNR > 8, we draw its kilonova light curve and calculate the value of P kn through the observation introduced above.The sums represent the expectation of the number of kilonova that WFST can detect through followup observations of BNS mergers.Table 1, 2 and 3 show the sum of P kn which normalized to one year's observation, with different EOSs, telescope's bands and exposure times.Here we list the results with the ms1 and wff2 models, which are the stiffest and softest EOSs mentioned in this paper.In addition, we pick the wff2 model as a comparison.As for observation bands, since the light curves for kilonovae in the u band change very rapidly and are less likely to be searched, WFST has a poor observation capability in the z band, we only consider single-band observations in the g, r, and i bands, and multi-band observations in the g and r/i bands in this work.Since the mass of the ejecta is related to the EOS of the NS, the stiffer EOS model ms1 corresponds to a lower observation probability, which are about 9.9 and 14.3 per year for the LHV and LHVKI networks, respectively, if the exposure time is 30s in the multi-band observations.The results are approximately the same for the bsk21 and wff2 models.
In the left panel of Figure 2, we compare the difference between the g band and i band P kn in the case of bsk21 model, 30s exposure time and the LHV network.For some of the kilonovae, the detection probability in the g band is significantly lower than in the i band, due to the light curves change more rapidly in the g band.Finally, the overall observation rates in the g band appear to be ∼ 20% and ∼ 30% lower than that in the i band for the LHV and LHVKI arrays, respectively.
Generally, a longer exposure times allow the telescope to have a larger observation depth, but on the other hand this will also decrease the area of the sky covered by the telescope when searching for kilonovae.Therefore, we also consider the results with two different exposure times, 60s and 90s.For these three cases, their readout time between each exposure is all set to 10 s.In the left panel of Figure 3, we compare the probability of each kilonova being observed with 30s and 90s exposure times.It can be seen that for most BNS mergers, increasing the telescope's single exposure time to 90s can help the search for kilonovae.In Table 1, 2 and 3, we compare the expectation of detection rates when using different exposure times.For the LHV network, using exposure times of 60s and 90s in multi-band observations will improve detection rates by ∼ 8% and ∼ 10%, respectively.As for the LHVKI array, it can be seen from Figure 1 that the improvement in the localization capability of the GW network will be significantly greater than the d L constraints.Therefore, the use of long exposure times facilitates the telescope to cover more of the observation volume.At this time, 60s and 90s exposure times improve the observation rate by ∼ 15% and ∼ 20%, respectively, compared to 30s exposure time.

NSBH mergers
Compared with BNS, NSBH will produce more powerful GW signals when it merges because of its larger mass.Therefore, for z < 0.2 NSBH mergers, about 30.5% and 50% of them can be detected by LHV and LHVKI network, respectively.We show the distributions of these observable NSBH mergers' redshifts, ∆Ω and ∆ ln(d L ) in the right panels of Figure 1.Since most of the observable NSBH mergers are at z > 0.1, the typical localization area given by LHV and LHVKI at 90% CL are ∼ 60 deg 2 and ∼ 20 deg 2 , which are a bit larger than the results of BNS mergers.The normalized detection numbers of LHV and LHVKI per year are 55 and 90, respectively.
According to Eq. ( 15), the spin of the BH will have a significant effect on the ISCO radius of the BH and further affect the matter ejection during NSBH merging.When χ BH deviates from 0, NSBH mergers will likely eject more matter and generate brighter kilonova emission.In general, the spins of binary mergers are related to the angular momentum transport in their stellar progenitors (Belczynski et al. 2020;Fuller et al. 2019;Fuller & Ma 2019), the environment where binary formed, and the tides and mass transfer during the orbit (Gerosa et al. 2018;Qin et al. 2019;Bavera et al. 2020).The results of GWTC-3 suggested that the effective spins of BH χ eff distribute around a very small value (The LIGO Scientific Collaboration et al. 2021b).However, if the binaries are located in the AGN accretion disk, they can gain spins by the accretion of surrounding gas (Yang et al. 2019;Safarzadeh et al. 2020;Hoy et al. 2022) and the observations of BH X-ray binaries suggest the existence of high spin BHs (McClintock et al. 2011;Miller et al. 2011;Qin et al. 2019).Thus, we follow the Zhu et al. (2021) and discuss the spin distribution for two special cases, χ BH ∼ N(0, 0.15 2 ) and χ BH ∼ N(0.85, 0.15 2 ) for the low-spin case and high-spin case, respectively.It is worth to mention that in the previous calculation we used the waveform for the nospinning system.Although the spin of the BH affects the waveform, it has no significant effect on the localization of the GW source (see Abbott et al. 2017c for example), so for the high-spin NSBH mergers, we used the same error bar as in the low-spin case for convenience.
In the low-spin case, we find that NSBH mergers can hardly produce strong kilonova emission.This situation is particularly evident in the soft EOS, for wff2 model, no NSBH merger can produce observable kilonovae in our simulations.For the stiffer EOS ms1, the multi-messenger observations rate for NSBH mergers is only ∼ 1/15 compared with BNS mergers.
And in the high-spin case, the NSBH mergers have greater probabilities to produce stronger kilonova emissions due to a smaller ISCO radius.For the stiffer model ms1, about 7 and 11 NSBH mergers per year can be multi-messenger observed for the LHV and LHVKI networks, respectively.Since the merging rate of NSBH is only about 1/3 that of BNS, these results are approximately 30% worse than BNS.In addition to χ BH , the EOS of the NS also has a large impact on the detection rate.In the case of bsk21 model, the observation rates are ∼ 4 and ∼ 7.2 per year for two networks.For the wff2 model, the detection rate of NBSH mergers is only ∼ 1/5 as much as the detection rate of BNS.
Since the kilonova emission from NSBH mergers is mainly contributed by the 'red' component, the light curves change slowly in all bands.Although the luminosity in g band is lower, WFST has a deeper detection depth in this band, so the detection probability of the kilonova emission from NSBH mergers in the g band and i band is approximately the same as seen in the right panel of the Figure 2.
For exposure time, due to a larger distribution for NSBH mergers' ∆Ω, it can be seen from the Figure 3 that, although the 90s exposure time still has a higher detection rate on the whole, more NSBH mergers have p 30s > p 90s compared with BNS systems.This suggests that for NSBH mergers, a shorter exposure time of 30s needs to be more favourably considered.
4.4.The applications in H 0 constraints Assuming that the other parameters of the ΛCDM model are determined, according to the propagation of uncertainty formula (Wall & Jenkins 2003), the measurement error of the Hubble constant can be written in the following form, where ∆ log d L comes from the measurement uncertainty of the GW detector networks, ∆z mainly comes from the peculiar velocity of the host galaxy, v p .In this work, we adopt σ v p ∼ 300 km s −1 .It is worthwhile to note that observations of kilonovae can place constraints on the inclination angle (Bulla 2019;Dhawan et al. 2020;Zhao et al. 2023), and then improve then improve the measurement for d L .However, these constraints on ι are highly dependent on models of kilonovae.Therefore, in this work we ignore the improvement of the d L measurement by observations of kilonovae.For N multi-messenger events, their constraining ability on H 0 is where p represents the probability that the i-th event's kilonova counterpart being observed.Through Eq. ( 25), we obtain the total constraints on the Hubble constant for these 10000 BNS/NSBH merger samples.When the statistical uncertainties are dominant, the uncertainty of H 0 is proportional to ∝ 1/ √ N, where N is the number of the GW events.After adopting the local merging rates of R BNSmergers,0 = 157.2Gpc −3 yr −1 and R NSBHmergers,0 = 48.7 Gpc −3 yr −1 and the redshift distribution model from Yu et   al. ( 2021), we predict that there will be 3000 and 900 z < 0.2 BNS, NSBH mergers, respectively, in a five years' observation time.The results normalized to these event numbers are listed in Table 4, 5, and 6.For BNS mergers, if the GW detector array is LHV and the kilonovae are searched in the r band for an exposure time of 90s, the H 0 can be constrained to ∆H 0 ∼ 2.8 km s −1 Mpc −1 .If the kilonovae need to be identified by multi-band observations, then the constraints become ∆H 0 ∼ 3.0 km s −1 Mpc −1 .In the case of the LHVKI array, due to the improved localization capability and d L constraints, these two results become ∆H 0 ∼ 1.7 km s −1 Mpc −1 and ∼ 1.9 km s −1 Mpc −1 , respectively.Therefore, multimessenger observations of BNS mergers will provide an important observational method for solving the Hubble tension.Chen et al. (2018) also predicted the ability of the 2G detector arrays to constrain the Hubble constant.They predicted that the LHV array could constrain the H 0 to ∼ 2% with two years of BNS multi-messenger observations, and ∼ 1% by adding two years of LHVKI's observations.After considering the effects of local merging rate, observation duration and duty cycle (0.5 for the LHV array and 0.3 for LHVKI in Chen et al. 2018), our constraints are still about three times looser.There are two main factors that account for this.First, Chen et al. (2018) assumed that all EM counterparts of BNS mergers detected by GW detectors would be observed.However, we find in our simulations that due to the influence of various factors, such as the large GW localization region, the Sun, the Moon, the Milky Waty, and the geographical location of the telescope, WFST can only be able to observe about 20%-25% kilonovae through the follow-up search.Secondly, Chen et al. (2018) predicted that the H 0 uncertainties would scale roughly as 15%/

√
N and 13%/ √ N for the LHV and LHVKI networks, respectively, where N is the detected BNS mergers number.In our calculation, these two results are 32%/ √ N and 26%/

√
N. This is due to Chen et al. ( 2018) chose 400 Mpc as the BNS detection threshold.In our calculations, we simulate the BNS mergers with z < 0.2 and find the LHVKI network can detect BNS mergers up to z ∼ 0.18.Considering that the redshift distribution of BNS merger is larger in our simulation, therefore the scaled H 0 uncertainties in our simulation are about two times larger.After taking these two factors into account, our results are roughly consistent with Chen et al. (2018).
For NSBH mergers, it is difficult to observe the kilonova emission in the case of low BH spin, so no effective constraints can be applied to H 0 .In the case of high spin, the results are significantly influenced by the NS EOS.For example, for the stiffer EOS ms1, the constraints on H 0 for the LHV and LHVKI are ∆H 0 ∼ 4 km s −1 Mpc −1 and ∼ 2.1 km s −1 Mpc −1 over a five -year observation period.The results with the LHVKI array are roughly comparable to those of BNS mergers.For the bsk21 model, the constraints become worse for ∆H 0 ∼ 5.4 km s −1 Mpc −1 and ∼ 2.8 km s −1 Mpc −1 in the cases of LHV and LHVKI arrays, respectively.And for the wff2 model, these two constarints are only ∆H 0 ∼ 8.2 km s −1 Mpc −1 for LHV and ∼ 4 km s −1 Mpc −1 for LHVKI.The results for bsk21 and wff2 models are significantly worse than the corresponding results for BNS.It shows that the contribution of multi-messenger observations of NSBH mergers to the Hubble constant constraints can reach the level of BNS mergers only if the BH spins are all very large and the EOS of the NSs are very stiff.However, this also means that multimessenger observations of NSBH mergers can impose strong constraints on BH spins and NS EOS.Therefore, the search for kilonovae produced by NSBH mergers will remain an important scientific goal for WFST in the future.Feeney et al. (2021) also investigated future multimessenger observations for NSBH mergers.With the A+ upgrade (The LIGO Scientific Collaboration et al. 2020), they predicted that about 2500 NSBH events will be detected by 2030.With the assumption of the DD2 NS EOS model, 99 of them will have sufficient ejecta (> 0.01 M ⊙ ) to produce observable EM counterparts.Feeney et al. (2021) predicted that multi-messenger observations of these NSBH mergers could reach 1.5-2.4% H 0 constraints.Considering that the mass of the NSBH merged ejecta is greatly influenced by the BH spin, mass ratio and the NS EOS, in some cases our results (low spin and/or soft EOS) differ dramatically from those in Feeney et al. (2021).Under the hard EOS model and highspin assumptions, our results are roughly similar.Our results together with those of Feeney et al. (2021) highlight the importance of improving our knowledge about the NSBH population.
It is worthwhile to note that for the BHs in NSBH mergers, we use a power law + peak mass distribution following the GWTC-3 results (The LIGO Scientific Collaboration et al. 2021b), which is mainly obtained from the BBH merging events.However, some recent studies, such as Broekgaarden et al. (2021) show that the BHs in NSBH systems may have a different mass distribution from those in BBH systems.Uncertainties in the mass distribution may introduce additional errors in the constraints of the Hubble constant.On the other hand, an accurate binary mass distribution model can be used to constrain cosmological parameters(The LIGO Scientific Collaboration et al. 2021b).Thus, It is important to obtain an accurate binary population model through future observations.We leave this to a future work.

4.5.
Uncertainties from fitting formulas In this section, we use several fitting formulas to calculate the mass of ejecta when BNS/NSBH merged.Since these formulas are fitted from tens of numerical relatively results with different NS EOS and merger parameters, and lack the constraints of the real observations, the final results depend on those fitting parameters and may have rather large uncertainties.In order to investigate the effect of the fitting parameters on the results, in this subsection we adopt several different fitting formulas for comparison.
For BNS mergers, we adopt the m dyn formula from Dietrich & Ujevic (2017), where fitting parameters are a = −1.
with a 1 = 0.007116, a 2 = 0.001436, a 4 = −0.02762,n 1 = 0.8636, and n 2 = 1.6840, and the numerical error is ∆M dyn = (0.1M dyn ) 2 + (0.01M ⊙ ) 2 .After replacing these fitting formulas, we repeat the above steps.For BNS mergers, most of the difference in ejecta mass is around 10%, with a fraction of mergers differing more in fitting results due to a lager mass ratio or NS mass.For NSBH mergers, the mass of the ejecta obtained by different fitting formulas is much more different due to large numerical error.Then take the case of a 30s exposure time for example, we show the multi-messenger observation rates and H 0 constraints with these fitting formulas in Table 7 and 8, respectively.For BNS mergers, the change in the multi-messenger observation rate is quite small, approximately a few percent smaller.However, for NSBH mergers, the difference in results is substantial, especially for ms1 and wff2, the two stiffest and softest models in our discussions with observation rates half as small and twice as larger, respectively.The difference under different fitting formulas further demonstrates the need for further studies of the kilonovae generated by NSBH mergers.
Reflecting on the H 0 constraints, the results of the BNS mergers are also nearly the same (about ≲ 4% looser), but the constraints of NSBH mergers with ms1 models are ∼ 1/3 looser.For NSBH mergers with wff2 models, five years' multi-messenger observations with the LHV network can constrain ∆H 0 ∼ 6.8 km s −1 Mpc −1 , which are 20% tighter that the results in Tabel 6.As for the LHVKI array, the constraints are only ∼ 10% tighter, this is due to that the additional part of the kilonovae mainly located at higher redshifts and contributing less to the H 0 constraints.

DARK SIRENS
In the previous section, we discuss the applications of BNS/NSBH merger-kilonova sirens in cosmology in the era of 2G GW detections.However, in the actual observation, the search for kilonova is limited by many factors, such as the Sun, the Moon, the Milky Way, the weather, and the geographical location of the telescope.For a ground-based telescope, its average observable area each night of the year is about 50% of the whole sky, and this percentage needs to be discounted considering the influence of weather, Moon and other factors.In particular, when a telescope is located at high latitudes, its observing sky area and observing hours in summer will be very limited due to the long hours of daylight.As a result, a significant fraction of kilonovae cannot be observed, even if they are quite bright or the GW signal is well localized.For example, in the simulations of the previous section, we find that for BNS mergers, the detection rate of their kilonova counterparts is only ∼ 1/5 of that of the LHVKI array.
In contrast, the observation of GW is hardly limited by these factors.When the GW detector is in operation, it can monitor GW signals coming from almost all directions.Therefore, for a dark siren, it is much less influenced by other factors than a bright siren.For example, the method of comparing GW events' localization area and a survey data is only affected by the localization accuracy of the GW signals and the completeness of the survey catalog in the localization region.Since the survey data are accumulated from long time observations, they are not affected by the constraints of the observation area, weather and other factors as much as follow-up observations.
In our pervious work Yu et al. (2020), we estimated the performance of the SBBH merger dark sirens by comparing the Sloan Digital Sky Survey data release 7 (SDSS DR7) group catalog (Yang et al. 2007) in the localization area.We found that for the 2G array, LHVKI, the effect of constraining the Hubble constant for the large mass SBBH dark sirens is close to that of the BNS bright sirens.However, for small-mass mergers, such as BNS, NSBH, and part of BBH systems, these events are difficult to be used as dark sirens to effectively constrain the Hubble constant in the era of 2G GW detection because of the weak GW signals produced.Only in the 3G era, with the significant increase in positioning capability, the BNS and NSBH dark sirens can play an essential role in cosmology.In this section, we will discuss the BNS and NSBH dark sirens and their applications.

Group catalog
Our group catalog is based on the high-resolution N-body Uchuu simulation (Ishiyama et al. 2021).Uchuu evolves the distribution of 12800 3 dark matter particles in a box of a side length of 2.0 h −1 Gpc and particle mass of 3.27 × 10 8 h −1 M ⊙ .The redshift range is from z = 127 to z = 0.The cosmological parameters adopted by Uchuu are listed at the end of Section 1.
For simplicity, we only use the halo/subhalo catalog obtained by the Uchuu simulation using the Rockstar halo finder (Behroozi et al. 2013) for the redshift z = 0 snapshot.We populate each subhalo in the catalog with a galaxy whose stellar mass is assigned according to the stellar to subhalo mass relation obtained in Yang et al. (2012).Once all the subhalos in the catalog are populated with galaxies of given stellar masses, we then calculate their accumulative stellar mass function.This accumulative stellar mass function is compared to the accumulative z-band luminosity function of galaxies obtained by Yang et al. (2021) from the DESI Legacy imaging Surveys data release 9. Using the abundance matching between the stellar mass and luminosity, we assign each galaxy a z-band luminosity.By properly taking into account the average k-correction (Yang et al. 2021) and the redshift space distortion effect (e.g, Yang et al. 2004), we construct a mock galaxy and group catalog with z-band apparent magnitude limit z limit = 21 in a sphere with redshift z ≤ 0.5 for this study (see Gu et al. 2023 in preparation for more details).
In Figure 4 we show the stellar mass distribution of galaxies and halo mass distribution of groups in catalog.The average stellar mass of galaxies and the average halo mass of groups are 2.10 × 10 10 M ⊙ /h and 1.94 × 10 12 M ⊙ /h, respectively.
To ensure the completeness of the group catalog, as in our previous work Yu et al. (2020), we choose halo mass M halo > 10 12 M ⊙ /h as a cutoff.These groups account for roughly 1/3 of the total stellar mass in the catalog.Under this halo mass cutoff, groups in the range z ≤ 0.5 are all roughly uniformly distribution in the comoving volume, which demonstrates the completeness of the catalog.We will discuss the reasonableness of adopting this cutoff at the end of this section.

Localization volumes
In our pervious work Yu et al. (2020), for a SBBH merger, we derived its covariance matrix Cov[α, δ, log(d L )] from the Fisher matrix and drew an ellipsoid in the parameter space of (α, δ, log(d L )), and converted the constraints on the source's luminosity distance d L to redshift z under the ΛCDM model.Then we traversed group catalog to find the number of groups in it.We used this number, N in to quantify the ability of identifying the source's host galaxy group.N in = 1 means that the host galaxy group of this GW event can be identified directly without EM counterpart.
Though for the CE2ET network, most SBBH mergers with z < 0.1 have N in = 1 due to a small localization volume (Yu et al. 2020), these results significantly rely on the cosmological model when tansforming d L to z.In contrast, the NSs' tidal effects in merging time provide a method of measuring GW sources' redshifts that does not rely on specific cosmological models.In this section, we discuss BNS and NSBH mergers respectively.For each group in catalog, we put an assumed BNS/NSBH merger at its center.Assuming that the EOS of NS has been tightly constrained from the X-ray plateaus (Lasky et al. 2014), kilonovae (Nicholl et al. 2021), etc., we obtain the covariance matrix Cov[α, δ, z] from 10-parameters' Fisher matrix, (α, δ, ψ, ι, M, η, t c , φ c , log(d L ), z), paint the error ellipsoid in the parameter space of (α, β, z) and traverse group catalog to count its N in . In

The distributions of N in
In each case of EOS, we calculate N in for the assumed GW events and get a distribution of N in .In Table 9 and 10, we list the fractions of N in = 1, N in ≤ 2, N in ≤ 5 and N in ≤ 10 for each case.
For BNS mergers, among these EOSs, ms1 has the highest fraction of N in = 1 due to smaller localization volumes, which reaches 44.8% and the fractions of N in ≤ 10 is 92.1%.In the case of ms1b, 43.5% BNS mergers' host galaxy groups can be identified directly through the observations of tidal effect.The results of other EOSs are worse than theirs, where in the cases of ap4, qmf40, sly, wff2 and pQCD800, the fractions of N in = 1 are below 30%.Meanwhile, in all cases, the fractions of N in ≤ 10 are larger than 70%.
For NSBH mergers, the results are much worse.Only for three EOSs, the fractions of N in = 1 exceed 10%, and they are H4, ms1 and ms1b.In the case of ms1, the fractions of N in = 1 and N in ≤ 10 are 16.1% and 49.3%, respectively.However, for wff2, only 3% NSBH mergers have N in = 1 and 15.4% assumed events have N in ≤ 10.
As a comparison, we also list the fractions of N in with error ellipsoids from Cov[α, δ, log(d L )] in Table 11.Almost all BNS, NSBH and BBH mergers have N in ≤ 10.And the fractions of N in = 1 are 53.5%,69.7%, 93.2% for BNS, NSBH and BBH mergers, respectively.We compare their probability density function (PDF) of N in with the cases of error ellipsoids from Cov[α, δ, z] in Figure 7 and 8.The ability to identify the BNS merger's host galaxy group through tidal effects with ms1 model is similar to the method of luminosity distance measurements.Thus for low redshift BNS mergers, when the EM counterparts are not observable, tidal effect observations will be a reliable and model-independent method of identifying the sources' host galaxy groups.However, for NSBH mergers, this method is not useful due to weaker tidal effect.
In Figure 9, we compare the host groups' redshift and stellar mass distributions for those N in events.From the left panels, we can see that BNS mergers around z ∼ 0.07 have the largest directly identified rate through the method of Cov[α, δ, log(d L )].In the case of tidal effect method with ms1 model, the distributions are almost same.However, in the case of wff2, the maximum rate occurs at z ∼ 0.06 due to a worse localization ability.For NSBH mergers, because of stronger SNR, the maximum N in = 1 rate in the case of Cov[α, δ, log(d L )] occurs at z ∼ 0.09.For the cases of Cov[α, δ, z] with ms1 and wff2 models, redshifts corresponding to the maximum rate are reduced to z ∼ 0.06 and ∼ 0.04, respectively.

APPLICATIONS OF DARK SIRENS IN H 0 MEASUREMENT
By comparing the galaxy/group catalog with GW signal's localization region, we can get GW event's redshift distribu- tion from the redshift statistics of its probable hosts.This approach was proposed by Del Pozzo (2012) and Chen et al. (2018), and has been applied in Fishbach et al. (2019) (2022) and Gray et al. (2022).Generally, these analyses use the constraints on {α, δ, log(d L )} given by the GW detector network, take specific value of Hubble constant H 0 as priors, transform the luminosity distance measurements into redshift space, compare them with the galaxy/group catalog, and calculate the likelihood of H 0 .Then, through the Bayes' theorem, we can obtain the posterior probability distribution of H 0 .Since this method relies on the prior of H 0 , the existence of galaxies/groups with different redshifts in the localized area will have a significant impact on the posterior of H 0 .
Observations of tidal effects, on the other hand, provide an independent method of measuring the redshifts of GW events.Thus, the search for potential hosts through constraints on {α, δ, z} will not rely on the prior of specific cosmological models.We will compare these two methods of dark sirens' redshift measurement in this section.
6.1.Constraints from Cov[α, δ, log(d L )] Following the Bayes' theorem, the posterior of H 0 can be written as where the GW likelihood p(d GW |d L (z, H 0 ), Ω) could be obtained from Cov[α, δ, log(d L )].We ignore the uncertainty of groups' positions and assume their redshifts follow a Gaussian distribution N(z obs , σ z |z), where z obs is the group's observed redshift, σ z is uncertainty of z obs .In Figure 10, we show the distribution of the difference between the groups' comoving redshift z com and observed redshift z obs .The distribution of ∆z = z obs − z com follows a Gaussian distribution with σ = 0.0014, hence we set σ z = 0.0014.Thus, the EM  likelihood p(d EM |z, Ω) could be written as where (z obs,i , Ω i ) represent the redshift and position of i-th possible host groups, and w i is the weight of each group, which represent the probability that different groups host a GW events.This weight is related to many factors, such as the stellar mass M s , morphologies, and metallicities (Cao et al. 2018).In this section, we choose two different weight function, w i = 1 and w i ∝ M s , as a comparison.
For the group prior p 0 (z, Ω), we assume the groups are isotropic distributed on large scales, and p 0 (z, Ω) could be written as (Chen et al. 2018;Fishbach et al. 2019) where the summation symbol is applied to the whole catalog.There is another approach that assume the galaxies/groups are uniformly distributed in comoving volume V. Therefore, p 0 (z, Ω) could be obtain from (Soares-Santos et al. 2019;Gair et al. 2023) where χ(z) is the comoving distance and H(z) is Hubble parameter with respect to redshift.In our previous work Yu et al. (2020), we discussed both approaches with the SDSS DR7 group catalog and found when the catalog is complete and the number density of the groups in comoving space does not vary with redshift, the results obtained by the two methods are almost the same.In this work, we are mainly concerned with the dark sirens at z < 0.1, a redshift range much smaller than the imcomplete redshift of our group catalog.Therefore, in this paper we only consider the former method.Following Cao et al. (2018) and Fishbach et al. (2019), the normalization term β(H 0 ) in Eq. ( 30) could be written as where P GW det (d L (z, H 0 ), Ω) and P EM det (z, Ω) are defined as where H is the Heaviside step function.For d EM , we set its thresold as z h < 0.5.For d GW , we set its threshold as SNR> 12. Since we only consider low redshift (z < 0.1) BNS/NSBH mergers, and the detection limit of 3G GW detector network is much farther than that (Yu et al. 2021), thus for all H 0 ∈ [20, 140] km s −1 Mpc −1 , the values of β(H 0 ) are same.
(37) 6.2.Constraints from Cov[α, δ, z] Since the observations of tidal effect can directly give constraints on z, the PDF of source's redshift could be obtained from the Bayes' theorem where p 0 (z) is the redshift prior.This expression can be expanded as where p 0 (Ω) represents the probability of producing GW events at Ω and p 0 (z, Ω) ≡ p 0 (z)p 0 (Ω).For this term, we set it to be proportional to the group distribution in catalog, so it can be obtained from Eq. ( 32).Using Eq. ( 31), Eq. ( 39) can be rewritten as In order to compare the abilities of these two methods in the Hubble constant constraints, we sample 200 assumed BNS/NSBH mergers at the center of the low redshift (z < 0.1) galaxy groups, where their redshifts equal to the groups comoving redshift z com and d L equal to d L (z com , H 0 = 70 km s −1 Mpc −1 ).In this section, to be consistent with the mock group catalog, we adopt the ΛCDM model and use the fixed value of Ω m mentioned in Section 1. First, we consider the case of w i = 1.The posterior distributions of H 0 from these two methods in this case are showed in Figure 11 and Table 12.The grey and blue lines in Figure 11 represent the constraints on H 0 for a single event and for 200 events, respectively.
In the upper left panel of Figure 11, we show the constraints from the method of Cov(α, δ, log d L ) and BNS mergers.It can be seen that since this method relies on the prior of H 0 , the posterior distribution given by a single event will exhibit a significant peak at values different from H 0 = 70 km s −1 Mpc −1 when there are other groups in the local- ization region.Therefore, when there are only a few events, these peaks may cause the measurement of the Hubble constant to deviate from the true value.However, as the number of events increases, the posterior distribution will approach a normal distribution with a central value of 70 km s −1 Mpc −1 .As a result, the constraint from these 200 BNS mergers is H 0 = 70.09+0.15 −0.15 km s −1 Mpc −1 .For NSBH samples, the result is H 0 = 70.10+0.14 −0.12 km s −1 Mpc −1 , which is slightly better than the case of BNS mergers due to a smaller error ellipsoid.The posterior distribution in this case is presented in the upper right panel of Figure 11.Now, we turn to the method of Cov(α, δ, z).The results in the case of ms1 model are show in the middle panels of Figure 11.For BNS mergers, since the identification of the host groups does not depend on the H 0 , the posterior distribution of all events is centered around 70 km s −1 Mpc −1 .Thus, when the number of events is small, the method of tidal effect observations will give a more reliable measurement of the Hubble constant.Eventually, the sum of the posterior distributions is H 0 = 69.99+0.15  −0.14 km s −1 Mpc −1 .This constraint is about 3% better than the constraint from the method of Cov(α, δ, log d L ).In the middle right panel, since a part of the NSBH mergers with larger mass ratio have a poor constraint on redshift, some peaks deviate from 70 km s −1 Mpc −1 .However, on the other hand, the stronger SNRs of the NSBH mergers lead to a smaller ∆ log(d L ).Thus, in total, the 200 NSBH mergers achieve a constraint of 70.05 +0.13  −0.13 km s −1 Mpc −1 for the Hubble constant, which is ∼ 10% better than BNS mergers.Since NSBH's local merging rate is one-third that of BNS's, after normalized, this constraint is ∼ 1/3 worse than the BNS mergers.
As a comparison, we also consider the EOS with the worst result in Section 5.3, wff2.The results are showed in the bottom panels of Figure 11.For BNS mergers, 200 events can constrain the Hubble constant to H 0 = 70.01+0.15  −0.15 km s −1 Mpc −1 .This result is ∼ 6% worse than the case of ms1 model.For NSBH systems, the constraint is almost the same with the case of ms1 model, which is H 0 = 70.06+0.13  −0.13 km s −1 Mpc −1 .This difference is significantly smaller than the difference in localization volume obtained by the two models.This is due to the fact that the measurement error in the d L dominates the H 0 constraints.The third and fourth columns represent the results when w i ∝ M s and w i = 1 weights are used in the calculation of the posterior respectively.
Usually, the choice of the weighting factor w i is a very important item in the dark standard siren method.When the number of groups in the localization volume is large, w i will have a large influence on the Hubble constant measurement and will make the results rely on the GW source population assumptions.
To investigate the effect of the weighting factor, similar to our previous work Yu et al. (2020), we consider a simple case where w i is proportional to group's stellar mass M s .We resample 200 BNS/NSBH mergers and repeat above calculations with w i ∝ M s .In addition, for comparison, we calculate the results for w = 1 with the same samples.The results are showed in Figure 12 and Table 13.It can be seen that the choice of two different weights has almost no effect on the error bar and only slightly changes the peak of the BNS mergers' constraints.This result is quite different from that of Gray et al. (2020).The main reason for this is that they discussed the 2G detector array, whereas we mainly focus on the 3G arrays.There is a several-order-of-magnitude difference between the localization volumes and the number of potential host galaxies/groups for GW events detected by the 2G and 3G arrays (Yu et al. 2020).When there are many potential host galaxies/groups, an incorrect binary population may cause a large bias in the posterior of the H 0 .
However, in our work, since most of the mergers have a small number of probable host groups, the contribution of these potential groups to the posterior of H 0 appears as a stacking of one or several isolated peaks, with the peak around H 0 = 70 km s −1 Mpc −1 corresponding to the 'true' host group, and the values of w i affecting the relative height of these peaks.For those 'golden' dark sirens with N in = 1, the choice of w i naturally has no influence, and due to the accuracy of the CE2ET array, they can provide very narrow constraints around H 0 = 70 km s −1 Mpc −1 .For other cases, the constraints given by those 'golden' dark sirens are equivalent to adding an extremely narrow prior to them around H 0 = 70 km s −1 Mpc −1 .Therefore the peaks from 'fake' host groups, which correspond to different H 0 , would be strongly suppressed or even eliminated.Since w i only affects the relative heights between each peak, it is a very minor quantity compared to the 'prior' produce by other samples.Therefore, we find roughly the same results for w i ∝ M s as for w i = 1.

The Gaussian likelihood approximation
In this paper, we use the Fisher matrix method, which is based on a linear approximation around the maximum likelihood and could forecast parameter errors in GW detection very simply and quickly.Recently, the reliability of the Fisher matrix method are widely discussed (Vallisneri 2008;Rodriguez et al. 2013;Grover et al. 2014;Magee & Borhanian 2022;Iacovelli et al. 2022;Gardner et al. 2023).They found that at low SNRs, since the linear approximation no longer holds, the likelihood function deviated from the Gaussian distribution.Thus the results obtained by Fisher's method are largely biased.However, when the SNR is very large (SNR > 25), due to the central limit theorem, these works found that the Fisher matrix and the Markov chain Monte Carlo (MCMC) methods have nealy the same estimations.For GW events with z < 0.1 in this section, their SNRs are much larger than 25, which guarantees the reliability of the localization area estimations.Therefore, to simplify the discussion, we consider that the results obtained by the Fisher matrix are the same as those of the MCMC approach.
In addition, to verify the effect of the variations in the likelihood function on the results, we make a simple test in Figure 13.It can be seen that the results for two different d L likelihood functions are almost the same, which somewhat indicate the reliability of the results.

Halo mass cutoff
In this section, we take a 10 12 M ⊙ /h group halo mass cutoff to ensure the completeness of catalog.To simplify the discussion, we assume that BNS/NSBH merging events are generated only in these groups.However, if a merging event happens in a small group outside of the catalog, the host group will be mislocated by the dark siren method, which may bias the H 0 measurements.Therefore, the reliability of the results needs to be tested.
For comparison, we continue assume that the distribution of mergers is exactly proportional to the stellar mass and resample the BNS/NSBH mergers with a cutoff throughout all groups without a halo mass cutoff.This way, there is a fraction of mergers whose host groups are not in the catalog.For a 10 12 M ⊙ /h group halo mass cutoff, this fraction is about 30%.First, to test the effect of mislocalization on the results, we extend the localization region to ensure that all mergers have at least one potential host.Due to these false localized events, the H 0 measurements of 200 BNS tidal effect dark sirens become 70.65 +0.15  −0.15 and 70.79 +0.16 −0.16 km s −1 Mpc −1 , for ms1 and wff2 model, respectively.It can be seen that the bias has exceeded 4 confidence intervals CIs.The reason is that when the host group of the BNS mergers is not in the catalog, its redshift measures will be determined by one or a few nearby large groups.Since the high redshift region has larger comoving volume, those 'false' host groups are more likely to occur at higher redshift, which lead to significantly larger results.And for NSBH systems, the biases are also lager than 4CIs.
Therefore, it is important to correct the bias caused by incomplete catalog.One method is to take use of the localization accuracy of the 3G detector network to exclude a part of mergers whose host groups are outside the catalog.For BNS mergers in the ms1 model, choosing a 90% localization volume cutoff can correct the H 0 measurements to 70.12 +0.17 −0.15 km s −1 Mpc −1 .However, in the case of wff2 model and NSBH mergers, the bias is almost unchanged due to the poor constraints in the redshift direction.
Another way is to complete the missing groups in the catalog.To account for the effects of these missing groups, we where G and Ḡ represent the host group is, and is not in the catalog, respectively.The first term can be obtained from Eq. (31).As for the latter term, it represents the contribution of these missing groups to the EM likelihood.Generally, this term is associated with the large-scale structure of the universe and the population of GW events.For simplicity, we neglect these effects and assume that p(d EM | Ḡ, z, Ω) is uniformly distributed.
Since groups with M halo > 10 12 M ⊙ /h occupy about 70% stellar mass, we set f comp = 0.7.
For 200 BNS mergers with ms1/wff2 model, the results are 69.98 +0.18  −0.18 /70.01 +0.19 −0.19 km s −1 Mpc −1 .It can be seen that after considering the contribution of these missing groups, the bias is nearly completely eliminated.In addition, since a flatter component is added to the redshift PDF for each merger, the final constraints on H 0 become slightly looser.For NSBH mergers, the results are 70.20 +0.16  −0.16 /70.20 +0.17 −0.17 km s −1 Mpc −1 .The difference between the results and the fiducial value is still less than 1.5 CIs.
Since the completeness factor relates to the fraction of stellar mass contained in the catalog, for comparison, we relax the halo mass cutoff to 10 11.5 M ⊙ /h, which makes the new catalog contain about 90% of the stellar mass.So with f comp = 0.9, the constraints become 70.07 +0.16  −0.16 /70.08 +0.17 −0.17 km s −1 Mpc −1 and 70.21 +0.15  −0.14 /70.20 +0.15 −0.14 km s −1 Mpc −1 , for BNS and NSBH mergers, respectively.It can be seen that with more groups contained in the catalog, the missing groups contribute much less to the EM likelihood, and therefore the constraints are tighter compared with the case with a 10 12 M ⊙ /h halo mass cutoff.
In brief, although in the results of Table 12 and 13, we neglected the mergers produced by small groups, which may lead to a bias in the measurement of H 0 , this bias can be well corrected by simply including the contribution of this part of missing groups in the EM likelihood.After adding these corrections, the mergers from the missing groups only slightly increase the error bar.Overall, a halo mass cutoff will not cause the results to be bias.

CONCLUSION
The detection of GW signals opens up a brand new approach to study the universe.Since the luminosity distance of a compact binary merger can be measured directly from waveforms, when its redshift is measured by other methods, this event can be used as a standard siren for cosmological research.
Among the various methods of redshift measurement, one of the most straightforward methods is to search for the EM counterpart of a GW event and identify the host galaxy to obtain the redshift information.In this paper, we use the Possis package to simulate the light curves of kilonova counterparts of low redshift BNS/NSBH mergers, then simulate multi-messenger observations of these samples by the 2G GW array together with WFST.For BNS mergers, WFST is expected to observe about 11 kilonovae per year by follow-up observations of the triggered signal from the LHV array.Over a five-year observation duration, these BNS bright sirens can constrain the Hubble constant to ∆H 0 ∼ 3 km s −1 Mpc −1 .As for the LHVKI array, since it can detect more BNS mergers and has higher localization accuracy, the multi-messenger detection rate of WFST is increased to about 18 per year, and the five-year observation time constraints are ∆H 0 ∼ 2 km s −1 Mpc −1 .
For NSBH mergers, its ejecta mass is mainly related to the relative size of the tidal disrupted radius and the ISCO radius, and is therefore strongly influenced by the mass ratio, NS EOS, and χ BH .Under the assumption that BH has high spin (χ BH ∼ N[0.85, 0.15 2 ]) as well as the NS EOS is the stiff ms1 model, NSBH mergers will have larger tidal disrupted radius, so the multi-messenger observation rates can reach ∼ 7 and ∼ 13 per year for the LHV and LHVKI arrays, respectively.For a middle EOS bsk21, these two rates reduce to ∼ 5 and ∼ 8 per year.For the softest EOS mentioned in the paper, wff2, the multi-messenger detection rates are only ∼ 2 and ∼ 4 per year left.For these three EOSs, the five-year observation duration constraints are ∆H 0 ∼ 2.1, 2.8, 3.9 km s −1 Mpc −1 with LHV, and ∆H 0 ∼ 4.0, 5.4, 8.2 km s −1 Mpc −1 with LHVKI, respectively.However, for the case of low BH spin, NSBH mergers can hardly produce bright kilonova emissions, thus WFST can hardly detect their EM counterparts.
In actual observations, most of the EM counterparts of compact binary mergers cannot be searched for due to various factors.Therefore, for these dark sirens, other methods are needed to measure the redshifts.Currently, a common method is to match the 3D localization area [α, δ, log(d L )] with a survey catalog to identify the host galaxies/groups of the GW source.In this method, it is necessary to use a d L − z relation to convert the constraints on d L into redshift.In the second half of this paper, we analyze the capability of the 3G GW detector array to identify the BNS/NSBH mergers' host group under the ΛCDM model with fixed parameters.With the [α, δ, log(d L )] constraints obtained from the Fisher matrix method, about 54% of the BNS and 70% of the NSBH host groups with z < 0.2 can be identified directly.However this approach relies on a specific cosmological model when transform d L to z.Therefore, we then discuss the feasibility of using the tidal effect of NSs at merging time to measure the redshift.Using the localization accuracy of third generation gravitational wave detector arrays and the redshift constraints obtained by measuring tidal effects, about 24% to 45% of the host galaxy groups of BNS mergers with z < 0.2 can be directly identified.The difference in this fraction is caused by the EOS of the NS, where the stiffer EOS corresponds to a higher fraction.Although the results are not as good as those of the d L -constrained method, this still indicates the feasibility of the method in identifying the host groups of BNS mergers.However, for NSBH mergers, the tidal effect is much weaker, so this fraction is only about 3%-16%.
It is worthwhile to note that the method of d L constraints needs to rely on cosmological model.Therefore, its measured redshift cannot be directly used for the H 0 measurement, but needs to be given a prior for the Hubble constant and then constrained using the Bayesian method.In this process, the error ellipsoid will be greatly magnified in the distance direction and causes the H 0 posterior given by a single event to exhibit multiple peaks that deviate from the actual value.The tidal effect approach, on the other hand, can give redshift constraints without relying on cosmological models.Benefiting from this, the tidally BNS and NSBH dark siren can constrain the H 0 to 0.2% and 0.3% over five years of observations, respectively, which are very close to the results obtained from the former method.
Del Pozzo et al. ( 2017) and Chatterjee et al. (2021) predicted that ET and CE can constrain the Hubble constant to ∼ 8% and ∼ 2% by measuring 10 3 BNS mergers' tidal effect.In this work, we consider array of 3G GW detector rather than a single interferometer.Due to the strong localization capability of the array, the accuracy of the redshift measurements can be significantly enhanced by comparing the localization volume with the galaxy group catalog.Thus our constraints reach ∼ 0.2% in 5 years' observation time, which are substantially better than those of Del Pozzo et al. (2017) and Chatterjee et al. (2021).In addition, for NSBH mergers, since their tidal effect is much weaker than that of BNS mergers, the direct measurement of redshift will be relatively poor.The method of comparing with catalog has greatly compensated this shortcoming and made the method of measuring the redshift of NSBH dark sirens by tidal effect become feasible.
An additional point to note is that in our simulations we use a mock catalog, which avoids the issues that exist in real surveys such as the edge effect.For the real catalog, if a part of the localized volume of the dark siren is outside the survey edge or the completeness range, those missing galaxies/groups may bias the H 0 measurement, so an integration of this part needs to be done (see Fishbach et al. (2019) for details).In particular, when the majority of the localization area is all outside the catalog range, the dark siren is not able to make an effective constraint on H 0 at this time.However, during the same period as the 3G GW detections (2030s), several spectroscopic experiments have been proposed to progress to Stage-5 spectroscopy for cosmology experiments, such as the MegaMapper (Schlegel et al. 2019), the Mauna Kea Spectroscopic Explorer (MSE; Marshall et al. 2019), and SpecTel (Ellis & Dawson 2019).We expect these experiments will provide a nearly complete group catalog, which will cover the majority of stellar mass at low redshift (Schlegel et al. 2022).

Foucart
et al. (2018) presented a fit of the remnant baryon mass outside of the BH after a NSBH merger based on 75 simulations,

Fig
Fig. 1.-The distributions of the redshifts, ∆Ω with 90% CL and ∆ ln(d L ) for BNS and NSBH mergers with SNR> 8.The left and right panels show the distributions of BNS and NSBH mergers, respectively.In each panel, the left y-axis represents the number distributions of 10000 samples obtained from the full populations of BNS/NSBH mergers, and the right y-axis N represents the distributions normalized to one-year's observation time.

Fig. 2 .
Fig. 2.-Comparison of the probabilities of observing kilonova between the g band and i band with the bsk21 model and 30s exposure time.The left and right panels show the distributions of BNS and NSBH mergers, respectively.

Fig. 4 .
Fig. 4.-The blue line represents the stellar mass distribution of galaxies in catalog and the red line represents the halo mass distribution of groups.
Figure 5, we show the cumulative distribution function (CDF) of localization volumes.We choose the cases of ms1 and wff2 as representatives, which have the smallest and largest localization volumes.As a comparison, we also show the CDF of localization volumes from the covariance matrix Cov[α, δ, log(d L )].For BNS mergers, the typical localization volumes from Cov[α, δ, log(d L )] are around (10 − 10 3 ) Mpc 3 .In the case of Cov[α, δ, z] and ms1, the localization volumes are slightly larger.And for wff2, the typical localization volumes are approximately 10 times larger.In the left panel of Figure 6, we show the distributions of BNS mergers' redshifts and their localization volumes obtained by different locating methods.At low redshift (z ≲ 0.05), the volumes from Cov[α, δ, log(d L )] are obviously smaller than those from Cov[α, δ, z].As the redshift increases, the tidal deformation term Φ tidal ( f ) will become larger.Thus, the localization capability of the observations of tidal effect will gradually catch up with the method of luminosity distance measurements.For NSBH mergers, the typical localization volumes from Cov[α, δ, log(d L )] are a few times smaller than BNS mergers, since they would produce stronger GW signals.However, in NSBH systems, the tidal effect is significantly weaker than in BNS systems due to the larger mass of the BH and its lack of tidal deformation.Therefore, the localization volumes from Cov[α, δ, z] are much larger.The distributions of NSBH mergers' redshifts and their localization volumes are showed in the right panel of Figure 6.Due to the existence of substructure in mass distribution of BHs around 34 M ⊙ , the localization volumes from Cov[α, δ, z] are divided into two parts.NSBH mergers with massive BHs will be located in a volume several magnitudes larger because of a weaker tidal field.

Fig. 5 .
Fig. 5.-The CDF of localization volumes for with BNS samples (blue line) and NSBH samples (red line).The soild line represents the case of Cov[α, δ, log(d L )].The dash and dot lines represent the cases of Cov[α, δ, z] with ms1 and wff2 models, respectively.
Fig. 6.-The distributions of localization volumes for different cases.The left-hand panel is for the case of BNS mergers and the right-hand panel is for the case of NSBH mergers.The red and blue dots represent the cases of Cov[α, δ, log(d L )] and Cov[α, δ, z] with ms1 model, respectively.

Fig
Fig. 7.-The PDF of BNS samples' N in from 1 to 20 with the CE2ET network.The left-hand panel shows the distribution with error ellipsoids from Cov[α, δ, z] and ms1 model.The right-hand panel shows the distribution with error ellipsoids from Cov[α, δ, log(d L )].
where d GW and d EM represent the GW signals observed by GW detector network and the data from galaxy group catalog, respectively.The term p(H 0 ) represents the prior of H 0 .In this section, we discuss a prior with uniform distribution in the interval[20, 140]  km s −1 Mpc −1 .Following Chen et al. (2018), the likelihood p(d GW , d EM |H 0 ) can be written as an integral over the solid angle Ω and redshift z p Fig. 9.-Redshift and stellar mass distributions of BNS mergers' host groups with N in = 1.The left and right panels represent the cases of BNS and NSBH mergers, respectively.The upper panels show the constraints from Cov(α, δ, log(d L )).The middle and lower panels show the constraints from Cov(α, δ, z).
40) where p(d GW |z, Ω i ) could be obtained from Cov[α, δ, z].As for the posterior distribution P(H 0 |d GW , d EM ), its CDF isP(H < H 0 |d GW , d EM ) = dz d L (z,H 0 ) 0 d dL p(z, dL |d GW , d EM ),(41) after differentiation of H 0 and combining the ∆ log(d L ) from the Fisher matrix method, we can obtain the posterior distribution of H 0 from P(H 0 |d GW , d EM ) = P(z|d GW , d EM )P(d L |d GW ) constraints of H 0 from different methods and EOSs with w i = 1.The numbers outside and inside the brackets denote the error bars for 200 mergers and normalized to five years of observation time, respectively.
Fig. 11.-The posterior distributions of Hubble constant H 0 from 200 BNS/NSBH samples with z < 0.1.The left and right panels represent the cases of BNS and NSBH mergers, respectively.The upper panels show the constraints from Cov(α, δ, log(d L )).The middle and lower panels show the constraints from Cov(α, δ, z).The grey lines show the constraints with a single event on H 0 and the blue line is the sum of them.
Fig. 12.-The same with Figure 11, but with w i ∝ M s .

Fig. 13
Fig. 13.-The comparison of the results with two different likelihood functions.The blue, red lines represent the results from the methods of Cov(α, δ, log(d L )) and Cov(α, δ, z), respectively.The solid and dashed lines reprsent the results of Gaussian likelihood functions about ln(d L ) and d L , respectively.All the constraints are obtained with ms1 model, w i ∝ m s , and 200 BNS mergers.modify the EM likelihood term p(d EM |z, Ω) by adding a completeness factor f comp (Gray et al. 2020), Multi-messenger observation rates of WFST at different bands for kilonovae per year when the NS EOS is wff2 model.

TABLE 9 The
fractions of the BNS samples with different values of N in and EOS.

TABLE 11 The
fractions of N in with error ellipsoids from Cov[α, δ, log(d L )].

TABLE 13 The
same with Table12, but resampled BNS/NSBH mergers with w i ∝ M s .