I N T RO D U C T I O N

The advent of gravitational wave measurements opened a new era of multimessenger 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 anticoincidence shield of the gamma-ray spectrometer on-board INTErnational Gamma-Ray Astrophysics Laboratory (Savchenko et al. 2017 ) within 0.1-0.647 s, and later complementing with the identification of the optical kilonova (e.g. Arcavi et al. 2017 ; E-mail: vivianeapa@cbpf.br Chornock et al. 2017 ;Coulter et al. 2017 ;Cowperthwaite et al. 2017 ;Evans et al. 2017 ;Kasliwal et al. 2017 ;Nicholl et al. 2017 ;Pian et al. 2017 ;Smartt et al. 2017 ;Soares-Santos et al. 2017 ;Tanvir et al. 2017 ;Valenti et al. 2017 ) detected about 11 h after the merger. This event produced the first direct and independent measure of H 0 , H 0 = 70 + 12 -8 km s -1 Mpc -1 (Abbott et al. 2017d ). After a 3-yr hiatus during which impro v ements 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 an ∼85 per cent impro v ement in sk y localization e xpected by the end of the run if all detectors are operating at their target sensitivities (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 an H 0 measurement with precision in the same order as what is achieved with other cosmological probes such as the cosmological microwave background (CMB; Planck Collaboration 2020 ) and the Cepheid (Riess et al. 2021 ) or Red Giant Branch (Freedman et al. 2019 )-MNRAS 528, 3249-3259 (2024) 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 (Verde, Treu & Riess 2019 ;Di Valentino et al. 2021 ). Despite these impro v ements, detection of the events electromagnetic counterparts remains a challenge, requiring dedicated follow-up 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 was the result of the coalescence of a 23.2 M sun black hole with a compact object of mass 2.5-2.67 M sun . Since the secondary mass lies in the mass mass-gap region, this object was either the massive neutron star or the lightest black hole ever seen in a binary system. Due to its excellent sky localization (23 degfoot_1 ), this event becomes a great candidate to provide the first detection of the counterpart of a binary system involving at least one black hole, hoping to shed light on the nature of this compact system. Several electromagnetic followups, from gamma rays to radio, were started by different groups (e.g. 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 ;Kilpatrick et al. 2021 ;Tucker et al. 2021 ;de Wet et al. 2021 ) with a continuous duration of up to more than 250 d after the merger. The properties of the electromagnetic counterpart candidates were analysed 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 gammaray 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 gamma-ray 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 H 0z de generac y and infer cosmological parameters. This methodology is known as dark standard sirens (see Gair et al. 2023 for a re vie w 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 H 0 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 Surv e y (DES) galaxy catalogue for binary black hole (BBH) mergers (GW170814 and GW190814, respectively) and showed that a single dark siren BBH provides a measure of H 0 with a precision of 48 per cent for GW170814 and 55 per cent 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 per cent against 18 per cent; Abbott et al. 2017d ).

Chen, Fishbach & Holz ( 2018 ) predicted that 5 years of detections for LIGO, Virgo and KAGRA collaboration (at design sensitivity) could lead to a precision of ∼ 5 per cent and 10 per cent of H 0 measurement for the BNS and BBH, respectiv ely. F or this result, the y assumed that all events within 10 000 Mpc 3 will be detected and that complete galaxy catalogues will be available. In the next decade, the arri v al 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 H 0 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 H 0 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 o v er, and because their localization re gion is co v ered by DELVEfoot_0 (Drlica-Wagner et al. 2021 ) galaxy catalogues. All the photometry redshift information is provided by the second release of DELVE data (DELVE DR2 2 , 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 photoz. This technique uses magnitudes and colour information to train the various Gaussian distributions that will be combined into the final PDF. In contrast to previous w ork, our w ork innovates by applying DELVE data to the standard siren methodology for the first time for two ne w e vents from the third observing run (O3), implementing a more refined artificial neural network technique for photoz measurements instead of the commonly used random forest algorithms (Zhou et al. 2020 ;Mucesh et al. 2021 ). The results of this study may provide insight into the potential of dark sirens as a cosmological probe, computing the precision level of H 0 measurement that this methodology can achieve with realistic photometric uncertainty and sky co v erage.

This paper is organized as follows: in Section 2 we describe the data used in the dark sirens methodology that is discussed in Section 3 . Our results are presented and discussed in Section 4 , and our final conclusions are presented in Section 5 . Throughout the article, we adopt a flat CDM cosmology with m = 0.3 and H 0 values in the 20-140 km s -1 Mpc -1 range. When not otherwise stated, quoted error bars represent the 68 per cent CI.

DATA

LIGO and Virgo data: gravitational wave events

Here, we extend the eight-event catalogue 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 second-half of the O3 (see Table 1 ), having a 90 per cent credible volume of 0.0034 Gpc 3 and a Table 1. Luminosity distance, 90 per cent CI area, and volume of gra vitational wa v e ev ents and candidates used in this analysis. We also report the reference paper or GCN that reports the sky map used for each e vent. These e v ents hav e 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 ).

The galaxies photoz's: a deep learning algorithm for DELVE data

The DELVE is a project that combines public data The photometric redshifts (or photozs) 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 MNRAS 528, 3249-3259 (2024) 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 photoz PDF, given some input features. The input features are the griz magnitudes, and the gr , gi , gz , ri , rz , iz colours. In the next sections, we use this approach to compute the photozs of the possible galaxy host whenever the spectroscopic redshift is not available.

The MDN was implemented with the following structure: a LMU layer with 212 units; a 2-layer Multi-Layer Perceptron with 96 units each; a Dropout layer with 20 per cent rate; and finally a Mix-tureNormal layer that returns the outputs (the mean, standarddeviation and weights of the 20 Gaussian distributions). The LMU layer was implemented using the keras-lmu Voelker, Kaji ć & Eliasmith ( 2019a ) application; the inner Perceptron and Dropout layers, the standard DL framework and the MixtureNormal output layer were built within the tensorflow and tensorflowprobability libraries API 3 (Abadi, Agarwal & et al. 2015 ).The architecture of the network also incorporates a Legendre Memory Unit (LMU, Voelker, Kaji ć & Eliasmith 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 photoz PDF estimated by the MDN output layer (also used for photometric redshift regression in the S-PLUS Survey in Lima et al. 2022 ) with the well-performing LMU layer to estimate the photometric redshifts. More details can be found in Teixeira et. al. Following the work of Zuntz et al. ( 2021 ), the LMU layer is included to more efficiently assign galaxies to redshift bins, selecting rele v ant information from previous data while simultaneously remo ving e xpendable data. F or 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 netw ork w as trained to maximize the PDF peak value for the spectroscopic redshifts ( z spec ) of each galaxy. The spectroscopic information came from a crossmatch between DELVE DR2 and the data available in different large sk y surv e ys (Colless et al. 2001 ;Mortlock, Madgwick & Lahav 2001 ;Wilson et al. 2006 ;Jones et al. 2009 ;Bacon et al. 2010 ;Drinkwater et al. 2010 ;Holwerda, Blyth & Baker 2011 ;Cooper et al. 2012 ;Mao et al. 2012 ;McLure et al. 2012 ;Bradshaw et al. 2013 ;F èvre et al. 2013 ;Newman et al. 2013 ;Baldry et al. 2014 ;Treu et al. 2015 ;Wirth et al. 2015 ;Bayliss et al. 2016 ;Momche v a et al. 2016 ;Nanayakkara et al. 2016 ;Tasca et al. 2017 ;McLure et al. 2018 ;Scodeggio, M. et al. 2018 ;Masters et al. 2019 ;Ahumada et al. 2020 ;Newman et al. 2020 ;Pharo et al. 2020 ;Mao et al. 2021 ;Mercurio et al. 2021 ), which resulted in approximately 4.5 million galaxies with z spec measurements. We also added the z spec 's available from the DECals DR9 Catalogue (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 arcsec.

In order to guarantee high quality photometric data used to train and test the model, we apply the following constraints on the colours, signal-to-noise ratio (SNR), and the limit of z spec :

(i) SNR > 3 for g 3 Tensorflow v2.9.1; Tensorflow Probability v0.17.0; keras-lmu v0.5.0

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 colour 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 colour-colour diagram in the regions -0.5 ≤ ri ≤ 1.5 and -0.5 ≤ ri ≤ 0.8. We restricted our z spec interval to a v oid spurious detections of low surface brightness galaxies located at high redshift. We also used the MODEST CLASS criteria (Drlica-Wagner et al. 2018 ) to remo v e contaminant stars by choosing the objects that lie in the classes 1 (high-probably 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 co v erage, with optical data in griz bands and partial co v erage when we are missing a band-co v erage only in (2) gri bands or (3) grz bands. Each MDN was trained with the magnitude (and colours) appropriate to the different scenarios. We used all of them to predict the z phot 's. The objects with full co v erage were assigned to the flag model GRIZ , and the same was made for gri and grz co v erages with the flags model GRI and model GRZ , respectiv ely. F or e xample, on the GW200202 154313 event we have approximately 3.4M objects with estimated z phot 's, being ∼ 62 per cent , ∼ 30 per cent and ∼ 8 per cent of the objects co v ered 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 31 252 and 6367 of these galaxies in the 90 per cent probability region of GW190924 021846 and GW200202 154313, respectively. Fig. 2 shows this final distribution, where we plotted the redshift distribution d N /d z subtracted from uniform number density (d N /d z) com assuming H 0 = 70 km s -1 Mpc -1 to emphasize the presence of o v erdensities along the line of sight.

To e v aluate the performance of our MDN method, we performed a complete analysis, e v aluating the point statistics and PDF's metric for the validation sample (represented by the 2.3 ×10 5 z spec 's that have not been used for training the photo-z's). The predicted photoz'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 + z spec )) results indicate that these data points only represent a minimum fraction ( < 4 per cent) of the entire sample o v er the redshift range of interest. In order to a v oid any systematic biases in DELVE galaxy distribution and their photoz, we select three different areas with the same size of LIGO 90 per cent probability region and analyse the photoz quality for these regions (solid lines and shadows in the right panel in Fig. 3 ). The right panel of Fig. 3 shows the median photoz bias in photoz 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 photoz bias is under control at z phot = 0.5, having median bias values smaller than 0.01 for each photoz bins and when considering the complete sample, the value reduces to -0.001. Thus, the measurements are uniform o v er the DELVE footprint. In contrast, the photoz 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 co v erage 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 + z spec , 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 < z spec < 0.3 and the σ 68 is less than 0.04 for all the photoz bins. These results are in agreement with previous each GW event, we start by computing the maximum redshift after converting the higher 90 per cent CI bounds in luminosity distance into the redshift, adopting the largest value of H 0 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 abo v e this threshold ( -19.39 and -20.32 for GW200202 154313 and GW190924 021846, respectively). The galaxies that survive this cut represent the most luminous galaxies. In this work we assume that the GW hosts trace the large-scale structure similarly to these most luminous galaxies. Thus, if the GW event occurred in a galaxy below our magnitude cut, the fainter galaxies follow the same matter distribution as the most luminous galaxies. The formation channel of binary systems has an influence on the specific properties of the GW hosts and whether it is able to pass through this cut. There are studies (e.g. Rauf et al. 2023 ) in the literature that point to the fact that more massive galaxies are indeed more likely to host BBH mergers, in fa v our of the chosen cuts. Ho we ver, more work on this front is needed to expand the discussion for the different formation channels and their parameter space (e.g. the assumptions that go in the binary population synthesis).

M E T H O D

In this work, we used the Bayesian formalism, described in detail in Chen, Fishbach & Holz ( 2018 ) and adapted into Soares-Santos et al. ( 2019) and Palmese et al. ( 2020Palmese et al. ( , 2023 ) ), to estimate the posterior probability of H 0 for the dark siren method. The H 0 posterior for a gra vitational wa ve measurement d GW and electromagnetic data d EM for a galaxy surv e y is written via Bayes' theorem as

where p ( H 0 ) is the prior on H 0 and p ( d GW , d EM | H 0 ) is the joint GW-EM likelihood. Assuming that the GW and EM measurements are independent, the joint likelihood can be written as p

. By marginalizing o v er the true redshift z, the sky position of the GW source, the photoz bias z and over all the possible galaxy hosts, the H 0 posterior can be written as in Palmese et al. ( 2023 ):

where r ( z, H 0 ) is the comoving distance,

is the Hubble parameter in a Flat CDM model, p ( z) is the prior on the photometric redshift bias which is measured with the method described in the Section 2.2 (see Fig. 3 ), β( H 0 ) is the selection function responsible for normalizing the likelihood, and Z is the evidence term defined as

is the marginal GW likelihood computed at the solid angle ˆ i and the redshift of the observed galaxy i , where we assume that this follow a Gaussian function according to Singer et al. ( 2016a ). The second term in the integral represents the marginal EM likelihood of the galaxy shifted by the photoz bias z, which is written as a product of galaxies photoz PDFs computed using the deep learning algorithm described in Section 2.2 . The abo v e expression includes the assumption that the source of the GW is located in one of the galaxies present in the galaxy catalogue, making it a function of the solid angle ˆ i and the redshift of each galaxy.

The abo v e posterior has two important ingredients: the selection effect defined by the β function and the photoz bias z. The first is associated with the selection effects adopted in the measurement process (of the electromagnetic counterparts and the detection of gra vitational wa ves), the β( H 0 ) function is computed following the same steps described in Chen, Fishbach & Holz ( 2018 ) and Palmese et al. ( 2023 ). For the electromagnetic emission selection effects, we used galaxies from the DELVE DR2 catalogue 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, Fishbach & Holz ( 2018 ) despite this being a simplification of the real EM selection effect, since it disre gards an y sk y accessibility, weather, and observing conditions, it is still a coherent approximation for estimating the observation of the real-time electromagnetic followup. On a large scale, we assume that galaxies are isotropically distributed across the sky. By marginalizing o v er the entire sky, the selection function can be written as

where p ( z) is the distribution of possible host galaxies and p GW sel ( d L ( z, H 0 ) ) is the probability of a source located at d L being detected. This term quantifies the GW selection effect introduced by detector sensitivity and detection conditions. For the computation of β( H 0 ) we follow the same steps as Palmese et al. ( 2023): first we simulate 70 000 BBH mergers for 20 different values of H 0 within our prior range [ 20 , 140 ] km s -1 Mpc -1 . 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 & Dickinson 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 & Price 2016 ;Singer et al. 2016a , b ) using the frequency domain approximant IMRPhenomD. Finally, we assume the O3 sensitivity curves for LIGO and Virgofoot_4 , use a matched-filter analysis, and calculate the SNR of each event. We assume, as a detection condition, that the network SNR is abo v e 12 and at least 2 detectors have a single-detector SNR above 4.

Another important effect considered in our analysis is the photoz 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 photoz, 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 photoz bias computationfoot_5 for the DELVE DR2 catalogue (see the detailed description of this calculation in Section 2.2 ) in different values of z that enter on H 0 posterior through p ( d EM | z, z).

This methodology can be extended to a sample of multiple events j with a combined data { d GW, j } and d EM , if we assume that the events are independent of each other and that they share the same galaxy catalogue. The Hubble constant posterior can be written as the product of the single event j likelihoods: Event

σ prior Reference GW190924 021846 70 . 4 + 54 . 7 -15 . 1 34.9 85 % This work GW200202 154313 51 . 2 + 61 . 6 -11 . 8 36.7 90 % This work GW170817-bright 68 . 8 + 17 . 30 -7 . 63 12.3 30 % Nicolaou et al. ( 2020 ) GW190814 78 + 57 -13 35 86 % Palmese et al. ( 2020 ) GW190924 + GW200202 60 . 33 + 55 . 79 -13 . 61 34.7 85 % This work GW190814 + GW170814 77 + 41 -22 31.5 77 % Palmese et al. ( 2020 ) 8 dark sirens 79 . 8 + 19 . 1 -12 . 4 15.8 39 % Palmese et al. ( 2023 ) 10 dark sirens 76 . 00 + 17 . 64 -13 . 45 15.55 38 % This work

Note that the abo v e e xpression does not apply to the two dark sirens studied here, since they do not share the same catalogue of galaxies because they are located in distinct regions of the sky and distance.

R E S U LT S A N D D I S C U S S I O N

We now use the DELVE photoz's in the dark siren methodology to produce the H 0 posterior for GW190924 021846 and GW200202 154313. Then we combine the results for these two ne w GW e vents with those for eight dark siren e vents (GW170608, GW170818, GW190412, S191204r, S200129m and S200311bg, GW170814, and GW190814) from Palmese et al. ( 2020Palmese et al. ( , 2023 ) ). The first five events were found in Palmese et al. ( 2023) using the DESI Le gac y Surv e y galaxies' redshifts, and the last two are presented in Palmese et al. ( 2020 ) with the photoz catalogue from DES. Fig. 4 shows the H 0 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 eight well-localized events (dark grey curve) found in Palmese et al. ( 2023 ). The two new events reduce the 68 per cent CI of the H 0 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 per cent and GW200202 154313 achieve the constraint of 90 per cent. The H 0 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 H 0 MNRAS 528, 3249-3259 (2024) (near 70 km s -1 Mpc -1 for 021846 and ∼51 km s -1 Mpc -1 for GW200202 154313) that is a consequence of the notable o v erdensity 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 o v er 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 H 0 ∼ 114 km s -1 Mpc -1 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 per cent) than GW200202 154313. The individual posteriors shown in Fig. 5 present a high probability at the high H 0 end. The same result was observed for the H 0 posterior of GW190814 and GW170814 in Palmese et al. ( 2020 ), as explained by the authors, this is a characteristic of the dark siren method once the GW analysis only provides a d L that is coherent with high values of H 0 . Here we decided not to adopt a wide H 0 prior, allowing to find a vast amount of galaxies distributed in redshifts that correspond to high values of H 0 . Fig. 2 reveals that the galaxies distributions for the two dark sirens are more uniformly distributed in comoving volume at high z, which results in a less informative EM likelihood and implies the reco v ery of the flat prior. 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 behaviour.

The combined result of all 10 dark sirens is shown in black in Fig. 4 . The mode of the final posterior and the 68 per cent CI is H 0 = 76 . 00 + 17 . 64 -13 . 45 km s -1 Mpc -1 . The addition of the two new dark sirens causes a reduction of ∼1 per cent o v er the 68 per cent confidence region found in Palmese et al. ( 2023 ). Table 2 summarizes our results and compares the performance with other standard siren analyses.

Fig. 6 illustrates the photoz bias effect on the H 0 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) H 0 values. For the GW190924 021846 dark siren, these v alues v ary 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 o v er the photo-z bias is minimum for all values of H 0 .

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 relati ve dif ference v arying between 0.0001 and ∼0.2 for the entire interval of H 0 .

In order to understand the impact of dark sirens on the precision of H 0 , we combine our results with the bright siren GW170817 from Nicolaou et al. ( 2020 ). Fig. 5 presents the combined H 0 posterior. The combination of GW190924 021846, GW200202 154313, and GW170817 gives H 0 = 68 . 84 + 15 . 51 -7 . 74 km s -1 Mpc -1 . We also combine the 10 dark sirens with GW170817, it gives H 0 = 71 . 54 + 10 . 96 6 . 61 km s -1 Mpc -1 representing an impro v ement of 6 per cent in the precision of the GW170817 measurement. This result highlights the impro v ement obtained when well-localized GW events at redshifts well co v ered by galaxy catalogues are incorporated into the analysis. Our results are in agreement with the recently presented results in Abbott et al. ( 2023 ), H 0 = 68 + 8 -6 km s -1 Mpc -1 , that used 47 dark sirens (43 BBH, 2 BNS, and 2 NSBH) from the third LIGO-Virgo-KAGRA GW transient catalogue (Collaboration et al. 2021 ) with GLADE + galaxy catalogues (D álya et al. 2018 , 2022 ). This constraint represents a reduction of ∼7 per cent in the 68 per cent CI of the H 0 measurement found with only GW170817. Although our measurements 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 H 0 . A possible cause of this dependency is the lack of completeness of GLADE + catalogue at the redshifts of interest.

C O N C L U S I O N S

In this work, we investigate the dark siren method to constrain H 0 , and present a new measure of H 0 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 photoz'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 and 51 . 2 + 61 . 6 -11 . 8 km s -1 Mpc -1 for GW190924 021846 and GW200202 154313, respectively, which is consistent with previous measurements of H 0 . The combination of GW190924 021846 and GW200202 154313 together with GW170817 bright siren leads to H 0 = 68 . 84 + 15 . 51 -7 . 74 km s -1 Mpc -1 , i.e. the addition of the two dark sirens reduces the 68 per cent CI interval by ∼7 per cent, which is comparable to the ∼12 per cent 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 surv e ys performing widesk y co v erage.

In addition, we also present the Hubble constant using only the dark standard siren method. We combine the H 0 posteriors found here with the posteriors of the eigh well-localized dark siren events (GW170814, GW190814, GW170608, GW170818, GW190412, S191204r, S200129m, and S200311bg) presented by Palmese et al. ( 2020Palmese et al. ( , 2023 ) ). The H 0 measurement found is 76 . 00 + 17 . 64 -13 . 45 km s -1 Mpc -1 , which has a precision of 20 per cent and the 68 per cent CI interval is ∼38 per cent of the prior width. Our result indicates that a sample with ten well-localized dark sirens and a complete galaxy catalogue can provide a significant constraint on the Hubble constant that is equi v alent 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, o v er 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 surv e ys, like the forthcoming Vera C. Rubin Le gac y Surv e y 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 H 0 at the several percentage level (Del Pozzo 2012 ). In this regime, we highlight the need for a more robust analysis, which takes into account potential systematics neglected in the methodology adopted here. Likely, some of the most significant sources of systematics will be the galaxy catalogue selection effects, galaxy catalogue 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 impro v e the dark siren methodology in order to consider these corrections.

AC K N OW L E D G E M E N T S

For the analysis, we use the PYTHON programming language, along with the following package: ASTROPY (Robitaille et al. 2013 ;Price-Whelan et al. 2018 ), MATPLOTLIB (Hunter 2007 ), NUMPY (Harris et al. 2020 ), SCIPY (Virtanen et al. 2020 ), TENSORFLOW (Abadi et al. 2015 ), and LIGO.SKYMAP (Singer & Price 2016 ;Singer et al. 2016a , b ;Kasliwal et al. 2017 ).

GT and CB acknowledge the financial support from CNPq (402577/2022-1). CB acknowledges the financial support from CNPq (316072/2021-4), FAPERJ (grants 201.456/2022 and 210.330/2022), and the FINEP contract 01.22.0505.00 (ref. 1891/22). The authors made use of Sci-Mind servers machines developed by the CBPF AI LAB team and would like to thank P. Russano and M. Portes de Albu-MNRAS 528, 3249-3259 (2024) querque for all the in infrastructure matters. CEM-V is supported by the International Gemini Observatory, a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America. GEM acknowledges support from the University of Toronto Arts & Science Postdoctoral Fellowship program. Time-domain research by DJS is supported by NSF grants AST-1821987, 1813466, 1908972, 2108032, and This manuscript has been authored by Fermi Research Alliance, LLC, under contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Go v ernment retains and the publisher, by accepting the article for publication, acknowledges that the United States Go v ernment retains a non-e xclusiv e, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Go v ernment purposes.