A dark siren measurement of the Hubble constant using gravitational wave events from the first three LIGO/Virgo observing runs and DELVE

The current and next observation seasons will detect hundreds of gravitational waves (GWs) from compact binary systems coalescence at cosmological distances. When combined with independent electromagnetic measurements, the source redshift will be known, and we will be able to obtain precise measurements of the Hubble constant H 0 via the distance-redshift relation. However, most observed mergers are not expected to have electromagnetic counterparts, which prevents a direct redshift measurement. In this scenario, one of the possibilities is to use the dark sirens method that statistically marginalizes over all the potential host galaxies within the GW location volume to provide a probabilistic redshift to the source. Here we presented H 0 measurements using two new dark sirens compared to previous analyses using DECam data, GW190924_021846 and GW200202_154313. The photometric redshifts of the possible host galaxies of these two events are acquired from the DECam Local Volume Exploration Survey (DELVE) carried out on the Blanco telescope at Cerro Tololo in Chile. The combination of the H 0 posterior from GW190924_021846 and GW200202_154313 together with the bright siren GW170817 leads to H 0 = 68 . 84 +15 . 51 − 7 . 74 km / s / Mpc . Including these two dark sirens improves the 68% confidence interval (CI) by 7% over GW170817 alone. This demonstrates that the inclusion of well-localized dark sirens in such analysis improves the precision with which cosmological measurements can be made. Using a sample containing 10 well-localized dark sirens observed during the third LIGO/Virgo observation run, we determine a measurement of H 0 = 76 . 00 +17 . 64 − 13 . 45 km / s / Mpc .


INTRODUCTION
The advent of gravitational wave measurements opened a new era of multi-messenger observation, shedding light on the properties of our universe.Standard sirens, a term introduced by Schutz (1986), provide a way to measure cosmological parameters by restricting the distance-redshift relation.The gravitational wave detections provide a direct measure of luminosity distance without any additional distance calibrator, justifying the name "standard sirens" in analogy with standard candles.If a source has an electromagnetic counterpart, its redshift (z ) can be directly measured, and we referred to them as "bright standard sirens".The first bright standard siren measured was the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017b), whose electromagnetic gamma-ray burst counterpart was detected by the Fermi Gamma-ray Burst Monitor (Goldstein et al. 2017) and the anti-coincidence shield of the gamma-ray spectrometer onboard INTErnational Gamma-Ray Astrophysics Laboratory (Savchenko et al. 2017) within 0.1 to 0.647 seconds, and later complementing with the identification of the optical kilonova (e.g.Arcavi et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Soares-Santos et al. 2017;Chornock et al. 2017;Kasliwal et al. 2017;Nicholl et al. 2017;Evans et al. 2017;Pian et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Valenti et al. 2017) detected about 11 hours after the merger.This event produced the first direct and independent measure of H0, H0 = 70 +12 −8 km/s/Mpc (Abbott et al. 2017d).After a three year hiatus during which improvements in the sensitivity of the detectors were made, the upcoming fourth run of the LIGO, Virgo and KAGRA collaboration will be able to observe a larger fraction of the universe than previous observing runs and projected to detect an estimated ∼90 gravitational wave events per year with a ∼85% improvement in sky localization of the source (Abbott et al. 2018).With more interferometers in operation (like the Einstein Telescope Sathyaprakash et al. 2012, Cosmic Explorer Abbott et al. 2017a and the LISA space interferometer Amaro-Seoane et al. 2017), it is possible that in the next years more standard sirens will be identified which can lead to a H0 measurement with precision in the same order as what is achieved with other cosmological probes such as the cosmological microwave background (CMB) (Planck Collaboration et al. 2020) and the Cepheid (Riess et al. 2021) or Red Giant Branch (Freedman et al. 2019)-calibrated type Ia supernovae.This new independent measurement of the Hubble constant can enable a way to clarify the origin of the observed current 4 − 6σ tension (Di Valentino et al. 2021;Verde et al. 2019).Despite these improvements, detection of the events electromagnetic counterparts remains a challenge, requiring dedicated followup campaigns and strategies (Bom et al. 2023), particularly for those events involving black hole companions, which may have no electromagnetic signature emitted or be associated to a flare (Bom & Palmese 2023;Rodríguez-Ramírez et al. 2023).
A prime example of the challenge to localize and identify the electromagnetic emission of an event involving a black hole is the GW190814 event (Abbott et al. 2020b).GW190814 has an excellent sky localization (23 deg 2 ) and a high probability (>99%) of being a black hole-neutron star merger, and became an excellent candidate to provide the first detection of the counterpart of these types of systems.
Several electromagnetic follow-ups, from gamma rays to radio, were started by different groups (e.g.Kilpatrick et al. 2021;Tucker et al. 2021;Dobie et al. 2019;Gomez et al. 2019;Ackley et al. 2020;Andreoni et al. 2020;Vieira et al. 2020;Watson et al. 2020;Alexander et al. 2021;de Wet et al. 2021) with a continuous duration of up to more than 250 days after the merger.The properties of the electromagnetic counterpart candidates were analyzed and compared with the theoretical prediction for NSBH fusion, including optical spectra, variability of radio sources, their location, photometric evolution, and redshift of possible host galaxies.Despite immense dedicated effort, no sign of a gamma-ray burst or any optical counterpart has been identified, but allowed to discard some possible types of electromagnetic transients such as: kilonova with large ejecta mass M ≥ 0.1M⊙ (Ackley et al. 2020), "blue" kilonovae with M > 0.5M⊙ (Kilpatrick et al. 2021), an AT2017gfo-like kilonova (de Wet et al. 2021), short gammaray burst with viewing angles less than 17 • (Kilpatrick et al. 2021), and a short gamma-ray burst-like Gaussian jet with a particular configuration (Alexander et al. 2021).In view of these problems, an alternative to the lack of an electromagnetic counterpart is to use the redshifts of galaxies that are within the coalescence location volume to break the H0-z degeneracy and infer cosmological parameters.This methodology is known as dark standard sirens (see Gair et al. 2023 for a review of the method).
The dark standard sirens approach was applied to constrain the cosmology in several LIGO and Virgo detections.Fishbach et al. (2019) studied the event GW170807 and showed that the obtained precision of H0 is about 3 times worse than the "bright" siren method (Abbott et al. 2017d).Soares-Santos et al. (2019) and Palmese et al. (2020) investigated the method with the Dark Energy Survey (DES) galaxy catalog for binary black hole (BBH) mergers (GW170814 and GW190814, respectively) and showed that a single dark siren BBH provides a measure of H0 with a precision of 48% for GW170814 and 55% GW190814.Recently, Palmese et al. (2023) demonstrated that 8 dark sirens well localized in the sky are able to provide a measurement as accurate as that obtained with a single bright siren GW170817 (about 20% against 18%; Abbott et al. 2017d).Chen et al. (2018) predicted that 5 years of detections for LIGO, Virgo and KAGRA collaboration (at design sensitivity) could lead to a precision of ∼ 5% and 10 % of H0 measurement for the BNS and BBH, respectively.For this result, they assumed that all events within 10,000 Mpc 3 will be detected and that complete galaxy catalogs will be available.In the next decade, the arrival of the next generation of terrestrial interferometers, such as the Einstein Telescope and the Cosmic Explorer, could rapidly increase the number of detections, allowing us to check the predictions of the percentage level of the measure of H0 made by Muttoni et al. (2023).
The intent of this study is to investigate the ability of the dark siren events GW190924_021846 and GW200202_154313 to constrain the Hubble constant.We combine our results with that of 8 dark sirens present in Palmese et al. (2023) and perform the most precise H0 measurement with the better localized dark sirens.The choice for these events is justified due to the small localization volume, which decreases the number of potential host galaxies to be marginalized over, and because their localization region is covered by DELVE1 (Drlica-Wagner et al. 2021) galaxy catalogs.All the photometry redshift information is provided by the second release of DELVE data (DELVE DR22 , Drlica-Wagner et al. 2022), the galaxy photometric redshift was estimated using the Mixture Density Network (MDN, Bishop 1994), a machine learning technique that provides the probability density function (PDF) of the photo-z.This technique uses magnitudes and color information to train the various Gaussian distributions that will be combined into the final PDF.In contrast to previous work, our work innovates by applying DELVE data to the standard siren methodology for the first time for two new events from the third observing run (O3), implementing a more refined artificial neural network technique for photo-z measurements instead of the commonly used random forest algorithms (Mucesh et al. 2021;Zhou et al. 2020).The results of this study may provide insight into the potential of dark sirens as a cosmological probe, computing the precision level of H0 measurement that this methodology can achieve with realistic photometric uncertainty and sky coverage.
This paper is organized as follows: in §2 we describe the data used in the dark sirens methodology that is discussed in §3.Our results are presented and discussed in §4, and our final conclusions are presented in §5.Throughout the article, we adopt a flat ΛCDM cosmology with Ωm = 0.3 and H0 values in the 20 − 140 km s −1 Mpc −1 range.When not otherwise stated, quoted error bars represent the 68% CI.

LIGO and Virgo data: gravitational wave events
Here, we extend the 8-event catalog used in Palmese et al. (2023) by adding two new events: GW190924_021846 and GW200202_154313.In total, our sample includes the 10 best localized events in the sky detected during the third LIGO/Virgo observing period.For these two added events, we used the gravitational wave data from the maps publicly available by the LIGO and Virgo collaboration in Abbott et al. (2021a) and Collaboration et al. (2021).The right ascension (RA), declination (dec), and distance probability are given in HEALPIX pixels (Górski et al. 2005), where this probability is supposed to be Gaussian along each line of sight.GW200202_154313 is the result of the merger of two black holes of approximately 7 and 10 solar masses, this is one of the best three-dimensional localizations from the secondhalf of the O3 (see Table 1), having a 90% credible volume of 0.0034 Gpc3 and a 90% CI area of 167 deg 2 .The location of maximum probability is centered at RA = 146.25 deg and dec = 20.98 deg.Marginalizing over all other parameters, the estimate of the luminosity distance has a mean equal to 364.3 Mpc with a standard deviation of 90.2 Mpc.The second detection, GW190924_021846 is the result of the merger of two lightest black holes (Abbott et al. 2021a), with inferred masses equal to m1 = 8.9 +7.0 −2.0 M⊙, m2 = 5.0 +1.4 −1.9 M⊙ and a 90% credible volume of ∼0.02 Gpc 3 and a 90% CI area of 348 deg 2 .The GW190924_021846 has a maximum probability of being located at RA, dec =(134.561,2.687) deg.At the peak location, the luminosity distance mean is 479.4Mpc and the standard deviation is 151.7 Mpc.The area of the sky enclosing 90% CI is 348 deg 2 .Fig. 1 shows the 90% CI contours of all events used in this work with the area covered by DELVE.The two new events are fully covered by DELVE, except for a small tip on GW190924021846 containing approximately 2%.The photometric redshifts (or photo-z s) for the DELVE data were computed using the deep learning method called Mixture Density Network.In brief, the method is a combination of a deep neural network with the assumption that any distribution can be written as a mixture of distributions (chosen to be the normal distribution in its traditional form).The deep neural network is trained, given some input features, to select the best parameters of the multiple distributions that will be mixed into a single distribution.The output parameters used are the mean, standard deviation, and mixing coefficients, which are the probabilistic weights of each normal distribution.In this way, the MDN is capable to reproduce the galaxy photo-z PDF, given some input features.The input features are the griz magnitudes, and the g-r, g-i, g-z, r-i, r-z, i-z colors.In the next sections, we use this approach to compute the photo-z s of the possible galaxy host whenever the spectroscopic redshift is not available.

The
The MDN was implemented with the following structure: a LMU layer with 212 units; a 2-layer Multi-Layer Perceptron with 86 units each; a Dropout layer with 20% rate; and finally a MixtureNormal layer that returns the outputs (the mean, standard-deviation and weights of the 20 Gaussian distributions).The LMU layer was implemented using the keras-lmu Voelker et al. (2019a) application; the inner Perceptron and Dropout layers, the standard DL framework and the MixtureNormal output layer were built within the tensorflow and tensorflow-probability libraries API 3 (Abadi et al. 2015).The architecture of the network also incorporates a Legendre Memory Unit (LMU, Voelker et al. 2019b) Layer at the head of the network.This architecture was one of the networks submitted in the LSST-DESC Tomography Optimization Challenge (Zuntz et al. 2021), and it exhibited the best performance for the DELVE DR2 photoz 's regression task.We combined the photo-z PDF estimated by the MDN output layer (also used for photometric redshift regression in the S-PLUS Survey in Lima et al. 2022) with Table 1.Luminosity distance, 90% CI area, and volume of gravitational wave events and candidates used in this analysis.We also report the reference paper or GCN that reports the sky map used for each event.These events have estimated false alarm rates of fewer than 1 in 10 3 − 10 23 years.These candidates have all recently been confirmed as gravitational wave events in Collaboration et al. (2021).2021), the LMU layer is included to more efficiently assign galaxies to redshift bins, selecting relevant information from previous data while simultaneously removing expendable data.For the loss function, we chose the maximum likelihood, which was minimized with the Nadam Optimizer (Dozat 2016) and results in a learning rate of 0.0002.The network was trained to maximize the PDF peak value for the spectroscopic redshifts (zspec) of each galaxy.The spectroscopic information came from a crossmatch between DELVE DR2 and the data available in different large sky surveys (Ahumada et al. 2020;Colless et al. 2001;Jones et al. 2009;Scodeggio, M. et al. 2018;Baldry et al. 2014;Mortlock et al. 2001;Newman et al. 2020;Drinkwater et al. 2010;Newman et al. 2013;Momcheva et al. 2016;Fèvre et al. 2013;Mercurio et al. 2021;Cooper et al. 2012;Holwerda et al. 2011;Mao et al. 2021;Bayliss et al. 2016;Bradshaw et al. 2013;McLure et al. 2012;Masters et al. 2019;Mao et al. 2012;Bacon et al. 2010;McLure et al. 2018;Wilson et al. 2006;Treu et al. 2015;Pharo et al. 2020;Tasca et al. 2017;Wirth et al. 2015;Nanayakkara et al. 2016), which resulted in approximately 4.5 million galaxies with zspec measurements.We also added the zspec's available from the DECals DR9 Catalog (Dey et al. 2019) also by doing a crossmatch with the DELVE DR2 data.All the matches were made considering a maximal separation of 0.972 arcseconds.
In order to guarantee high quality photometric data used to train and test the model, we apply the following constraints on the colors, signal-to-noise ratio (SNR), and the limit of zspec: The SNR cuts were used to eliminate spurious sources, bad measurements, and very faint galaxies.The g mag limitation serves to reinforce the exclusion of faint galaxies.The color cuts were made in order to eliminate nonphysical (extremely blue and extremely red) objects (see Drlica-Wagner et al. 2018), thus the majority of the objects in our sample populate the color-color diagram in the regions −0.5 ⩽ r − i ⩽ 1.5 and −0.5 ⩽ r − i ⩽ 0.8.We restricted our zspec interval to avoid spurious detections of low surface brightness galaxies located at high redshift.We also used the MODEST_CLASS criteria (Drlica-Wagner et al. 2018) to remove contaminant stars by choosing the objects that lie in the classes 1 (highprobably galaxy) and 3 (ambiguous classification).
To account for the lack of a band on our data, we decided to train 3 different MDN's for each different observation scenario: (1) full coverage, with optical data in griz bands and partial coverage when we are missing a band -coverage only in (2) gri bands or (3) grz bands.Each MDN was trained with the magnitude (and colors) appropriate to the different scenarios.We used all of them to predict the z phot 's.The objects with full coverage were assigned to the flag model_GRIZ, and the same was made for gri and grz coverages with the flags model_GRI and model_GRZ, respectively.For example, on the GW200202_154313 event we have approximately 3.4M objects with estimated z phot 's, being ∼ 62%, ∼ 30% and ∼ 8% of the objects covered by griz, grz, gri bands, respectively.
After the selection cuts, our training sample contains about one million objects distributed at redshift z < 1.There are 31252 and 6367 of these galaxies in the 90% probability region of GW190924_021846 and GW200202_154313, respectively.Fig. 2 shows this final distribution, where we plotted the redshift distribution dN /dz subtracted from uniform number density (dN /dz )com assuming H0 = 70 km/s/Mpc to emphasize the presence of overdensities along the line of sight.
To evaluate the performance of our MDN method, we performed a complete analysis, evaluating the point statistics and PDF's metric for the validation sample (represented by the 2.3×10 5 zspec's that have not been used for training the photo-z's).The predicted photo-z 's as a function of the measured spectroscopic redshifts is shown in the left panel of Fig. 3.We can see that the majority of data points lie close to the diagonal, thus pointing to the accuracy of the predicted redshifts.Additionally, we can see the presence of outliers in every redshift interval.However, the outlier fraction (which is defined as |∆z| > 0.15 × (1 + zspec)) results indicate that these data points only represent a minimum fraction (< 4%) of the entire sample over the redshift range of interest.In order to avoid any systematic biases in DELVE galaxy distribution and their photo-z, we select three different areas with the same size of LIGO 90% probability region and analyze the photo-z quality for these regions (solid lines and shadows in the right panel in Fig. 3).The right panel of Fig. 3 shows the median photo-z bias in photo-z bins of size 0.025 for DELVE and LEGACY-DR9 measurements.The results for DELVE full spectroscopic sample (dashed red line) and DELVE limited areas (dashed red line) revealed that the photo-z bias is under control at z phot = 0.5, having median bias values smaller than 0.01 for each photo-z bins and when considering the complete sample, the value reduces to -0.001.
Thus, the measurements are uniform over the DELVE footprint.In contrast, the photo-z results from the LEGACY full spectroscopic sample (blue dashed line) appear to outperform DELVE, displaying median bias values consistently below 0.005.This difference in quality could be attributed to the fact that LEGACY measurements benefit from uniform coverage across all bands and also leverage the advantages of infrared bands in their Spectral Energy Distributions (SEDs).The scatter of z phot predictions was quantified with the normalized median absolute deviation, defined as σNMAD = 1.48 × median (|∆z|/ (1 + zspec)), and the 68th percentile width of the bias distribution about the median (σ68).Data for DELVE objects brighter than r < 21 yields σNMAD = 0.023 for all the galaxies in 0 < zspec < 0.3 and the σ68 is less than 0.04 for all the photo-z bins.These results are in agreement with previous works that use similar techniques to measure photometric redshift in large sky surveys (see, Lima et al. 2022).In order to validate the individual photo-z PDF as a whole, we use two different metrics: the probability integral transform (PIT) distribution and the Odds value.The PIT distribution for DELVE data has a positive skewness, which indicates that our deep learning methods overestimate the z phot .The Odds value represents the fraction of the photo-z PDF that is contained in the interval zspec ±0.06, its distribution reveals that our MDN models produce narrow photo-z PDFs with their values centered near 1.
As shown previously, (Soares-Santos et al. 2019;Palmese et al. 2020;Palmese et al. 2023) to overcome the fact that our sample is magnitude-limited, we have to ensure that it is volume-limited.For that, we follow the same steps used in Palmese et al. (2023).For each GW event, we start by computing the maximum redshift after converting the higher 90% CI bounds in luminosity distance into the redshift, adopting the largest value of H0 we considered in the prior.The next step is to find an absolute magnitude threshold value that corresponds to the apparent magnitude limit at the maximum redshift.Finally, we exclude all galaxies in our sample that have an absolute magnitude above this threshold (-19.39 and -20.32 for GW200202_154313 and GW190924_021846, respectively).

METHOD
In this work, we used the Bayesian formalism, described in detail in Chen et al. ( 2018) and adapted into (Soares-Santos et al. 2019; Palmese et al. 2020;Palmese et al. 2023), to estimate the posterior probability of H0 for the dark siren method.The H0 posterior for a gravitational wave measurement dGW and electromagnetic data dEM for a galaxy survey is written via Bayes' theorem as p (H0|dGW, dEM) ∝ p (dGW, dEM|H0) p (H0) , where p (H0) is the prior on H0 and p (dGW, dEM|H0) is the joint GW-EM likelihood.Assuming that the GW and EM measurements are independent, the joint likelihood can be written as p (dGW, dEM|H0) = p (dGW|H0) p (dEM|H0).By marginalizing over the true redshift, the sky position of the GW source, the photo-z bias ∆z and over all the possible galaxy hosts, the H0 posterior can be written as in Palmese where r (z, H0) is the comoving distance, H (z) = H0 Ωm (1 + z) 3 + 1 − Ωm 1/2 is the Hubble parameter in a Flat ΛCDM model, p (∆z) is the prior on the photometric redshift bias, β (H0) is the selection function responsible for normalizing the likelihood, and Z is the evidence term defined as Z = dzip (dEM|zi) r 2 (zi) /H (zi).The above expression includes the assumption that the source of the GW is located in one of the galaxies present in the galaxy catalog, making it a function of the solid angle Ωi and the redshift of each galaxy.
The above posterior has two important ingredients: the selection effect defined by the β function and the photo-z bias ∆z.The first is associated with the selection effects adopted in the measurement process (of the electromagnetic counterparts and the detection of gravitational waves), the β (H0) function is computed following the same steps described in Chen et al. ( 2018) and Palmese et al. (2023).For the electromagnetic emission selection effects, we used galaxies from the DELVE DR2 catalog distributed up to the known absolute magnitudes for each of the GW events in the analysis, where we consider only those in z < 0.5.As reported in Chen et al. ( 2018) despite this being a simplification of the real EM selection effect, since it disregards any sky accessibility, weather, and observing conditions, it is still a coherent approximation for estimating the observation of the real-time electromagnetic follow-up.On a large scale, we assume that galaxies are isotropically distributed across the sky.By marginalizing over the entire sky, the selection function can be written as where p (z) is the distribution of possible host galaxies and p GW sel (dL (z, H0)) is the probability of a source located at dL being detected.This term quantifies the GW selection effect introduced by detector sensitivity and detection conditions.For the computation of β (H0) we follow the same steps as Palmese et al. (2023): first we simulate 70,000 BBH mergers for 20 different values of H0 within our prior range [20,140] km/s/Mpc.The BBH population is distributed through the redshift distribution p (z) which is a function of the merger rate evolution and the cosmologydependent comoving volume element.For simplicity, we assume that the merger rate follow the Madau-Dickinson star formation rate (Madau 2014).The mass of the black holes is distributed according to a power-law with index 1.6 (in agreement with the results found in Abbott et al. 2021b).We draw spins from a uniform distribution between (−1, 1).The GW signals were generated using the BAYESTAR software (Singer 2016;Singer et al. 2016b;Singer et al. 2016a) using the frequency domain approximant IMRPhenomD.Finally, we assume the O3 sensitivity curves for LIGO and Virgo4 , use a matched-filter analysis, and calculate the SNR of each event.We assume, as a detection condition, that the network SNR is above 12 and at least 2 detectors have a single-detector SNR above 4.
Another important effect considered in our analysis is the photo-z bias correction.When we are dealing with simulated data, the machine learning algorithm used for photometric redshift estimates can provide a biased redshift probability distribution function.The non-uniform training samples can cause systematic biases in the photo-z, causing the peak of the distribution to be shifted by ∆z from the true value of z.In order to consider this effect, we use the photo-z bias computation5 for the DELVE DR2 catalog (see the detailed description of this calculation in section 2.2) in different values of z that enter on H0 posterior through p (dEM|z, ∆z).
This methodology can be extended to a sample of multiple events j with a combined data {d j GW , d j EM }, if we assume that the events are independent of each other and that they share the same galaxy catalog.The Hubble constant posterior can be written as the product of the single event j likelihoods: p (H0|{dGW,i}, dEM) ∝ p (H0) p (dEM|H0) j p (dGW,j|H0) (4)

RESULTS AND DISCUSSION
We now use the DELVE photo-z 's in the dark siren methodology to produce the H0 posterior for GW190924_021846 and GW200202_154313.Then we combine the results for these two new GW events with those for 8 dark siren events (GW170608, GW170818, GW190412, S191204r, S200129m and S200311bg, GW170814 and GW190814) from (Palmese et al. 2020;Palmese et al. 2023).The first five events were found in Palmese et al. (2023) using the DESI Legacy Survey galaxies' redshifts, and the last two are presented in Palmese et al. (2020) with the photo-z catalog from DES. Fig. 4 shows the H0 posterior from the combination of these two new dark sirens (dark red curve) and the final result (black curve) after combining the posterior of all the ten dark siren events.For comparison, we also show the results (blue curve) found in Palmese et al. (2020) with the dark sirens GW170814 and GW190814 and for the 8 well-localized events (dark gray curve) found in Palmese et al. (2023).The two new events reduce the 68% CI of the H0 prior to values close to those found in Palmese et al. (2023) (see Table 2): GW190924_021846 is able to reach the value of 85% and GW200202_154313 achieve the constraint of 90%.The H0 posterior distributions for GW190924_021846 and GW200202_154313 are presented in Fig. 5, we can see that both dark sirens display an evident peak at a low value of H0 (near 70 km/s/Mpc for GW190924_021846 and ∼51 km/s/Mpc for GW200202_154313) that is a consequence of the notable overdensity of galaxies (see Fig. 2) around redshift 0.05 to 0.1 and 0.05 for GW190924_021846 and GW200202_154313, respectively.As a result of better localization volume, which corresponds to a marginalization over a smaller number of galaxies, we can see that the posterior of GW200202_154313 has a narrower peak, but the presence of a secondary peak at H0 ∼ 114 km/s/Mpc makes it flatter than GW190924_021846 (its kurtosis value is lower, ∼1.79, than that produced by GW190924_021846, ∼1.95).The analysis of the skewness showed that event GW190924_021846 produces a slightly more asymmetric posterior (a relative difference of approximately 58%) than GW200202_154313.The individual posteriors shown in Fig. 5 present a high probability at the high H0 end.The same result was observed for the H0 posterior of GW190814 and GW170814 in Palmese et al. (2020), as explained by the authors, this is a consequence of the fact that the GW analysis only provides a dL estimate that can correspond to either a low z and H0 or a high z and H0.The cut on the prior range does not bias the final result, as it only changes the redshift range considered in the dark siren analysis and not the posterior behavior.
The combined result of all 10 dark sirens is shown in black in Fig. 4. The mode of the final posterior and the 68% CI is H0 = 76.00 +17.64  −13.45 km/s/Mpc.The addition of the two new dark sirens causes a reduction of ∼1% over the 68% confidence region found in Palmese et al. (2023).Table 2 summa-rizes our results and compares the performance with other standard siren analyses.
Fig. 6 illustrates the photo-z bias effect on the H0 posterior distribution for the events GW190924_021846 and GW200202_154313.We see that the effect is a little more significant for GW200202_154313, with a relative difference of ∼0.1(0.09) for low(high) H0 values.For the GW190924_021846 dark siren, these values vary from 0.0002 to 0.18.This is different than what was discussed in Palmese et al. (2020), which showed that the effect of marginalization over the photo-z bias is minimum for all values of H0.
Another correction applied here is the full redshift PDF instead of a Gaussian approximation.The effect of this correction is shown in Fig. 7, where it is almost the same for the two dark sirens, with the relative difference varying between 0.0001 and ∼ 0.2 for the entire interval of H0.
In order to understand the impact of dark sirens on the precision of H0, we combine our results with the bright siren GW170817 from Nicolaou et al. (2020).Fig. 5 presents the combined H0 posterior.The combination of GW190924_021846, GW200202_154313, and GW170817 gives H0 = 68.84+15.51  −7.74 km/s/Mpc, which is in agreement with the recently presented results in et al. ( 2023), H0 = 68 +8 −6 km/s/Mpc, that used 47 dark sirens (43 BBH, 2 BNS, and 2 NSBH) from the third LIGO-Virgo-KAGRA GW transient catalog (Collaboration et al. 2021) with GLADE+ galaxy catalogs (Dálya et al. 2018(Dálya et al. , 2022)).This constraint represents a reduction of ∼7% in the 68% CI of the H0 measurement found with only GW170817.Although our measurement is less precise than Abbott et al. (2023) (as expected given the smaller number statistics), we note that we expect our result to be less sensitive to the black hole population assumptions.As noted in Abbott et al. (2023), these assumptions, and specifically the shape of the mass distribution, strongly dominate the inference on H0.A possible cause of this dependency is the lack of completeness of GLADE+ catalog at the redshifts of interest.By combining these results with the 8 dark sirens from Palmese et al. (2023), we find an improvement of 6% in the precision of the GW170817 measurement.This result highlights the improvement obtained when well-localized GW events at redshifts well covered by galaxy catalogs are incorporated into the analysis.

CONCLUSIONS
In this work, we investigate the dark siren method to constrain H0, and present a new measure of H0 provided by two GW events detected by LIGO/Virgo, GW190924_021846 and GW200202_154313, with the redshifts of the potential host galaxies derived using DELVE DR2 data.The estimation of galaxies photo-z 's was performed using the deep learning technique Mixture Density Network.Our analyses implement the full redshift PDF of the galaxies instead of the Gaussian approximation.The main result of this study includes the measurement of the Hubble constant of 70.4 +54.7  −15.1 km/s/Mpc and 51.2 +61.6 −11.8 km/s/Mpc for GW190924_021846 and GW200202_154313, respectively, which is consistent with previous measurements of H0.The combination of GW190924_021846, GW200202_154313, together with GW170817 bright siren leads to H0 = 68.84+15.51  −7.74 km/s/Mpc, i.e. the addition of the two dark sirens reduces the 68% CI interval by ∼7%, which is comparable to the ∼12% found in Palmese et al. (2020) when they add GW190814 and GW170814.This result demonstrates the power of well-localized dark siren events in better constraining the determination of the Hubble constant using deep imaging photometry obtained from surveys performing widesky coverage.
In addition, we also present the Hubble constant using only the dark standard siren method.We combine the H0 posteriors found here with the posteriors of the 8 well-localized dark siren events (GW170814, GW190814, GW170608, GW170818, GW190412, S191204r, S200129m and S200311bg) presented by Palmese et al. (2020) and Palmese et al. (2023).The H0 measurement found is 76.00 +17.64  −13.45 km/s/Mpc, which has a precision of 20% and the 68% CI interval is ∼38% of the prior width.Our result indicates that a sample with ten well-localized dark sirens and a complete galaxy catalog can provide a significant constraint on the Hubble constant that is equivalent to that achieved with a standard siren, providing complementary information to the standard method.
Our results provide an indication of the dark siren potential as a precision cosmological probe.After a period of sensitivity upgrades, over the past few months, the LIGO/Virgo/KAGRA collaboration has returned to operation and is expected to make ∼ 90 detections of mergers per year (Abbott et al. 2018).With the increase in GW observations and the arrival of deeper and wider surveys, like the forthcoming Vera C. Rubin Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), it is possible that in the next few years, dark sirens will provide a measure of H0 at the several percentage level (Del Pozzo 2012).In this regime, we highlight the need for a more robust analysis, that takes into account potential systematics neglected in the methodology adopted here.Likely, some of the most significant sources of systematics will be the galaxy catalog selection effects, galaxy catalog completeness, the dependence of host galaxy properties on the BBH formation channels, and the use of a Gaussian approximation for the GW likelihood instead of its full asymmetric distribution.In future work, we intend to improve the dark siren methodology in order to consider these corrections.Table 2. Hubble constant measurements found with the dark sirens from the three LIGO/Virgo runs and the bright siren GW170817.All priors are flat in the range [20,140].The uncertainty from the flat prior is derived by assuming the same H 0 maximum found in the analysis.Quoted uncertainties represent 68% HDI around the maximum of the posterior.The "σ H 0 /σ prior " column shows the 68% CI from the posterior divided by 68% CI of the prior width.versity, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Am-

Figure 1 .
Figure 1.LIGO/Virgo GW dark standard sirens analyzed in this paper, where the contours represent the 90% CI localization from the sky maps.The shaded regions are those that are covered by the DELVE catalogs used in this work.
Abbott et al. (2021a) the well-performing LMU layer to estimate the photometric redshifts.More details can be found inTeixeira et.al.  (in preparation).Following the work ofZuntz et al. (

Figure 2 .
Figure2.Redshift distribution of galaxies in the 90% CI area of the dark siren events analyzed in this work.The distribution is subtracted with a dN/dz with uniform number density to highlight the presence of overdensities and underdensities along the line of sight.The dashed blue line shows the distribution using the photometric redshift point estimates from the DELVE, the red line shows the same redshifts when their uncertainty is considered as a Gaussian error.The gray vertical lines represent the luminosity distance of each GW event marginalized over the entire sky, assuming an H 0 of 70 km/s/Mpc, and the shaded regions are the 1σ uncertainties considering the same H 0 ; these regions are only shown for reference.

Figure 3 .
Figure 3. Photometric redshift quality assessment plots using the testing sample from the available spectroscopic data.Left: density plot of galaxies in the validation sample, showing the predicted z phot 's (PDF peaks) as a function of spectroscopic redshift.The red dotted lines represent the outliers limits, where outliers z out phot are defined as |z out phot − zspec| > 0.15/(1 + zspec).Right: median value of the bias distribution ∆z = z phot − zspec in bins of photometric redshift for our model used in the DELVE DR2 and the Legacy Survey DR9 photometric redshifts.The median was calculated taking 3 different non-correlated regions in the sky covered by our test sample.The regions were chosen to have the same area as GW200202_154313.Both panels are plotted with r < 21 mag, zspec < 0.3 and z phot < 0.3.

Figure 4 .
Figure 4. Hubble constant posterior distributions for the dark sirens considered in this work and previous works.The black line is the result of the combination of the two dark sirens considered here with the posteriors of all the 8 dark sirens (gray line) found in Palmese et al. (2023).The combination H 0 posteriors of GW190924_021846 and GW200202_154313 is shown in red, and for comparison we also show the combination result of GW170814 and GW190814 (blue line) presented in Palmese et al. (2020).For comparison, we show the 1σ constraints on H 0 found by Planck Collaboration et al. (2020), Riess et al. (2021) (R22) and Abbott et al. (2023) (GWTC-3) as the vertical shaded regions.

Figure 5 .
Figure 5. Hubble constant posterior distributions found using the DELVE galaxies for GW200202_154313 (green line) and GW190924_021846 (blue line).The dashed line represents the GW170817 bright siren result adapted from Nicolaou et al. (2020), which includes the peculiar velocity corrections for the galaxy host NGC 4993.The black line is the result from the combination of the two dark sirens from this work with GW170817, and the vertical dashed lines show the 68% region for this posterior.Posteriors are arbitrarily rescaled only for visualization purposes.

Figure 6 .
Figure 6.Photo-z bias effect on the Hubble constant posterior distributions for GW190924_021846 (blue lines) and GW200202_154313 (yellow lines).Solid curves are the results considering the inclusion of the photo-z bias correction, and dashed curves ignore this correction.The sub-plots present the relative difference between these two curves.

Figure 7 .
Figure 7.Comparison between the Hubble constant posterior distributions for GW190924_021846 (blue curves) and GW200202_154313 (yellow curves) obtained by using the full galaxies redshift PDFs and a Gaussian approximation.The subplots present the relative difference between the curves.
H 0 (km/s/Mpc) σ H 0 (km/s/Mpc) σ H 0 /σ prior and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State Uni- Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey.The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l'Espai (IEEC/CSIC), the Institut de Fisica d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.Based on observations at Cerro Tololo Inter-American Observatory, NSF's NOIRLab (NOIRLab Prop.ID 2019A-0305; PI: Alex Drlica-Wagner), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.BASS is a key project of the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program "The Emergence of Cosmological Structures" Grant # XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance.The BASS is also supported by the External Cooperation Program of Chinese Academy of Sciences (Grant # 114A11KYSB20160057), and Chinese National Natural Science Foundation (Grant # 11433005).This manuscript has been authored by Fermi Research Alliance, LLC, under contract No. DE-AC02-07CH11359 with the US Department of Energy, Office of Science, Office of High Energy Physics.The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.
paro, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnolo-gia e