A dark standard siren measurement of the Hubble constant following LIGO/Virgo/KAGRA O4a
C. R. Bom V. Alfradique A. Palmese G. Teixeira L. Santana-Silva A. Santos P. Darc








0000-0003-4383-2969
]C. R. Bom
Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil
Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Rodovia Márcio Covas, lote J2, quadra J - Itaguaí (Brazil)
Clecio R. Bom
debom@cbpf.br

0000-0002-5225-1923]V. Alfradique
Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil

0000-0002-6011-0530]A. Palmese
McWilliams Center for Cosmology, Department of Physics Carnegie Mellon University, 5000 Forbes Ave, Pittsburgh, 15213, PA, USA

0000-0002-1594-208X]G. Teixeira
Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil

0000-0003-3402-6164]L. Santana-Silva
Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil


Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil

0000-0001-8833-474X]P. Darc
Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil






We present a new constraint on the Hubble constant (H_0) from the standard dark siren method using a sample of  5 well-covered gravitational wave (GW) alerts reported during the first part of the fourth LIGO/Virgo/KAGRA observing runs in combination with standard dark sirens from the first three runs. Our methodology relies on the galaxy catalog method alone. We use the full probability density estimation of photometric redshifts derived by a deep learning method using the DESI Legacy Survey and DELVE galaxy catalogs. We add the constraints from the binary black hole mergers candidates S231226av, S231206cc, S230919bj, S230627c, and S230922g to the sample of standard dark sirens analyzed in <cit.>. We combine the H_0 posterior for 5 new standard sirens with other 10 previous events (3 with updated posteriors), finding H_0 =  69.9^+13.3_-12.0  km s^-1 Mpc^-1 (68% Highest Density Interval) with the catalog method alone.  This result represents an improvement of ∼ 23% comparing the new 15 dark siren constrain with the previous 10 dark siren constraint, and a reduction in  uncertainty of ∼ 40% from the combination of 15 dark and bright sirens compared with the GW170817 bright siren alone.
The combination of dark and bright siren GW170817 with recent jet constraints yields H_0 of 68.0^+4.3_-3.8  km s^-1 Mpc^-1,  a ∼ 6% precision from Standard Sirens, reducing the previous constraint uncertainty by ∼ 10% . 















§ INTRODUCTION


The current 4-6σ tension in the Hubble constant <cit.> arises from the significant discrepancy between different cosmological probes, in particular from the Cosmic Microwave Background <cit.>, and those using Supernovae (SN) and Cepheids for the local distance ladder <cit.>. New independent measurements of the Hubble constant (H_0) have the potential to shed light on this discrepancy (e.g., ), and, depending on their precision, could arbitrate the tension. Among novel probes, the standard sirens <cit.> methodology employs gravitational wave occurrences to obtain luminosity distances. This information is used to infer cosmological parameters, most notably the Hubble constant  H_0, upon integration with redshift data derived from host galaxies. This emergent probe could play an important role as it is independent of the cosmic distance ladder <cit.>.






Standard sirens are categorized as "bright" in instances where an electromagnetic counterpart is definitively identified alongside a singular host galaxy, and as "dark" or "statistical" in the absence of such counterparts. Bright sirens, exemplified by GW170817 <cit.>, yield measurements of high precision. However, the requirement of host identification poses a series of challenges due to the wide search volume and the cadence requirements <cit.>. Also, kilonovae have been the only widely confirmed electromagnetic sources detectable from gravitational waves, although black hole mergers are proposed to produce electromagnetic counterparts in certain circumstances with a few identified candidates <cit.>. Therefore, the dark siren method can be applied to a larger number of events, including GW170817 <cit.>, the binary black hole mergers GW170814 <cit.> and GW190814 <cit.>, and several events from the first three LIGO/Virgo observing runs <cit.>. Dark standard sirens assuming the catalog method <cit.> rely on the position and redshift of potential host galaxies, leading to less precise results than bright sirens on a single-event basis. However, combining dark and bright sirens can enhance constraints on cosmological parameters, leveraging the abundance of events without counterparts.






The fourth LIGO/Virgo/KAGRA (LVK) Observing run (O4) began on May 24, 2023, and is scheduled to span 20 months, including two months allocated for commissioning breaks dedicated to maintenance.
The initial segment of O4 referred to as O4a, finished on January 16, 2024. During O4a there have been 82 GW alerts reported, of which 80 were classified as binary black hole candidates. It is important to note that the Virgo detector was not  operational during the O4a run, therefore the typical sky localizations were worse than would have been expected if it had joined the run. Meanwhile, KAGRA has participated in the run for a limited duration, albeit at a notably reduced binary neutron star (BNS) inspiral range compared to the LIGO detectors.

In this work, we use a set of well-localized and confident candidate events from O4a covered by public photometric catalogs and imaging data, mainly by the Dark Energy Spectroscopic Instrument (DESI) Legacy Survey <cit.> and DECam Local Volume Exploration Survey <cit.> , derived the photometric redshifts using full probability density estimation of photo-z using the same approach from <cit.> to adding 5 new sirens expand the sample of dark sirens from <cit.> to 15 sirens. The dark siren catalog method used in this work to constrain H_0 is described in the aforementioned papers, we also refer to the Appendix <ref> for a description of our methodology.

























§ DATA









 §.§ The LIGO/Virgo GW data



In this study, we expand the analysis upon the previous works with 8-events <cit.> and 10 events <cit.> from LVK run O1-O3, using the catalog method. We employ comparable selection criteria from prior studies (see P23 and A24) and, as such, investigate five new sirens identified in O4 with  > 70% of their probability coverage falling within the Legacy or DELVE surveys, and with luminosity distances less than d_L<1500 Mpc. 

We restricted in distance based on the previous works, as the catalogs turn more incomplete and also the dependency on Ω_m becomes more relevant. The positional data for these events comes from the maps publicly provided by the LVK collaboration, we present the 90% Credible Interval (CI) sky region of all events used in this work in Fig. <ref>. These maps right ascension (RA), declination (dec), and distance probability, are represented using HEALPIX pixelation. Within this framework, the probability distribution along each line of sight from the maps is assumed to follow a Gaussian distribution. We draw attention to the fact that the criteria adopted here are chosen with the intention of selecting those superevents with greater constraint capability, but there is no impediment for others to be added to the sample, as long as the selection function (see its definition in Appendix <ref>), defined in the H_0 posterior, correctly describes the cuts adopted.





 The five novel sirens candidates events we selected are S231226av <cit.>, S231206cc <cit.>, S230922g <cit.>, S230919bj <cit.> and S230627c <cit.>. They are all classified as Binary Black Hole mergers (BBH) with probability > 99%, except for S230627c, which has a classification between neutron star-black hole (49%), BBH (48%) with (3%) of being noise. In particular, S231226av  is among the lowest False-Alarm-Rate (FAR) events, while S230919bj and S230627c are in the top 20% percentile of low FAR in O4. The novel events 90% volume are comparable to the previous sirens studied with 90% volume ∼ 10^-3 Gpc^3. We use the latest public skymaps (see Table <ref> for details) produced by the python code Bilby <cit.>, that are available on the GraceDB event page [<https://gracedb.ligo.org/>]. We also use the last skymaps from the GW Transient Catalog (GWTC [<https://gwosc.org/eventapi/html/GWTC/>], ) for three events from P23 which are updated on this work namely GW191204_171526, GW200129_065448 and GW200311_115853.






 §.§ Optical Survey Data


We made use of publicly available catalogs from the DESI Legacy Survey <cit.> and DELVE survey <cit.>. The combined footprint is presented in Fig. <ref>. We use this data to obtain precise photometric redshifts using the same technique and Deep Learning method described in A24. We use Mixture Density Network (MDN) to derive the full probability density functions for each galaxy in the survey catalog. A more comprehensive description of the model and photometric redshift quality assessment in both surveys and a comparison between the MDN method and the public data for the Legacy Survey is presented in Appendix <ref>. Our final constraints use the Legacy Survey Probability Density Function (PDF) instead of Gaussian approximations from the Legacy Survey photo-z catalogs <cit.>. The Legacy Survey results slightly outperformed the DELVE for the same model. We use DELVE-based photometric redshift as a validation and we found a small impact for the final constraints of < 0.5 km/s/Mpc.


























































































    





In the GW sample processed in this work, only the S230919bj event has significant 90% region  uncovered by the photo-z catalogs from both Legacy galaxy survey and DELVE. To compensate for the insufficient coverage, we adopted the same procedure from <cit.>. Our strategy involves distributing simulated galaxies in regions lacking data. 

To ensure that marginalization includes all possible host galaxies and leave our Hubble constant measurement free of an underestimated uncertainty, the injected galaxies follow our prior distribution that is given by the training sample. The photo–z distribution of these fake galaxies was sampled by the Monte Carlo technique. We assume a uniform spatial distribution, and the number density is given by the value of the Legacy Survey galaxy catalog. For the photo–z precision, we first found the relation between the photo-z error and photo-z by computing the mean and standard deviation value of the photo-z error in photo-z bins of size 0.05 from our training sample, after that, we sampled from a Gaussian function for each photo-z bins. This same procedure was adopted for the apparent magnitude distribution of these fake galaxies, since it is intrinsically associated with the redshift value. Here, we also assume that all the fake photometric redshift PDFs follow a Gaussian function; we make a Gaussian sampling around the assigned true value with the standard deviation given by their respective photo-z error. In section <ref> we present the impact of these fakes galaxies and found it to be minor.



§ RESULTS AND DISCUSSION
 



In this section, we present the H_0 posterior produced using the dark siren methodology described in the Appendix <ref> with the Legacy Survey photo-z’s PDFs constructed using the deep learning methodology from A24 and also described in <ref>. We computed the H_0 posterior for five GW candidate BBH events detected in the LIGO O4a run and reanalyzed with the updated skymaps tree events from P23, as they were not yet confirmed as gravitational wave events (see <ref>). The 68% CI of the five O4a superevents represent about 69% to 91%, depending on the event analyzed, of the 68% CI of the prior width.
The most constraining event from O4a is S230627c, because of its better localization, galaxy catalog coverage, and quality of the photo-z at the relevant redshift ranges. For the revisited gravitational wave events, despite the increase in 90% CI volume, all the three events previously studied with the  maps in P23 present an improvement of ∼ 18%,  5% and 9% in 68% CI, respectively.  

The main cause of these improvements is the quality of our photo-z measurement (see Fig.  <ref>), in Appendix <ref> we discuss these results in more detail. 

The majority of the H_0 posteriors show a clear peak, except for S231226av and GW191204, which exhibit multiple peaks in different ranges of H_0 (raging from 49 to 136 km/s/Mpc). As already mentioned in previous work () these peaks are associated with the presence of overdensities regions along the line of sight (see Fig.<ref> in Appendix <ref>). In the dark siren methodology, this is indicative of the host galaxy being more likely to be in these regions. Additionally, it is evident that certain events display a more prominent peak, which may stem from the presence of a prominent overdensity within an improved localization. In other words, the need for marginalization over fewer galaxies and/or the fact that the galaxies live at similar redshift, lead to a more informative H_0 posterior.


The H_0 posterior for GW191204  presents a considerable probability at H_0≈ 130 km/s/Mpc, which is associated with an overdensity of galaxies at z≈0.25 (see Fig.<ref> in Appendix <ref>). We also can see that this H_0 posterior begins to decrease at the high–H_0 end, due to the absence of galaxies at high redshift (the galaxy density in Fig.<ref> is negative at the high z end, indicating the presence of an underdense region).

The combination of the 15 dark sirens from O1-O4a is depicted in Fig.<ref> by the green line and detailed in table <ref>. The maximum a posteriori and the 68% CI is H_0 = 69.9^+13.3_-12.0 km/s/Mpc.
As it can be seen, our results are consistent within 1σ with Planck <cit.> and the local Cepheid-supernova distance ladder <cit.>. This result agrees with the latest dark sirens study  using the BBH population and catalog method <cit.>, that found H_0 = 67^+13_-12 km/s/Mpc for 47 dark sirens from the third LIGO-Virgo-KAGRA GW transient catalogue <cit.> and the  galaxy catalogue.  A possible justification for the similar uncertainties despite the difference in number of the events used is the completeness and inhomogeneity of the  catalog at the redshift of interest of the studied dark sirens. This contributes to the dependence between the final H_0 constrain and the BBH population assumptions observed in <cit.>, and minimizes the galaxy method contribution to the final constraint. Furthermore, we also combine our dark siren results with the bright siren GW170817 from <cit.>, which considers the effect of peculiar velocity to the H_0 constrain first found by <cit.>. From this analysis, we find H_0 = 69.20^+8.98_-5.85 km/s/Mpc. The dark siren information provides a notable reduction of 41% in the 68% CI and a 7% improvement in the relative precision. We further combine our final dark siren H_0 posterior with the bright siren analysis for GW170817 presented in <cit.>, which corrects the peculiar velocity contribution in the measure of H_0 presented in <cit.> where the superluminal motion measured by the Very Large Baseline Interferometry (VLBI, ) and the afterglow data is used to measure the inclination angle of the GW170817 event. The addition of an independent EM observation measure helps break the degeneracy between the distance and the inclination angle, which is one of the main contributions to the uncertainty of the gravitational wave distance measurement. We find that the dark sirens information cause a ≈ 12% improvement of the 68% CI, and leads to our final constraint: H_0=68.00^+4.28_-3.85 km/s/Mpc. Furthermore, we combine our 15 dark sirens result with a more recent H_0 measurement <cit.>, which combines the GW170817 measurements with the electromagnetic counterpart associated with afterglow rather than the superluminal jet motion as done in <cit.>. This result presents an agreement of 1.4σ with the H_0 measurements of Planck <cit.>, which indicates a reduction of 0.1σ of the original result.  


To guarantee the use of the most appropriate galaxies catalog in terms of coverage of the localization volumes of the gravitational wave events treated here and the quality of the estimated photo-z's, we also analyzed the measurements of H_0 for our dark sirens sample only using the DELVE catalog. This allows us to evaluate if the H_0 measurement and photo-z quality outperform our results obtained using the Legacy survey. 

The photometric redshifts for the DELVE galaxies were computed with the deep learning methods presented in A24 adapted for the redshift range of interest here. The results of the individual  H_0 posterior distribution found with the Legacy and DELVE galaxy catalog are in agreement, the combination of the 15 dark sirens found with DELVE leads to a higher H_0 peak value (approximately 2.0 km s^-1 Mpc^-1) and a lower uncertainty of H_0 (σ_H_0∼ 0.5 km s^-1 Mpc^-1) compared to those presented in table <ref>. This dark siren result combined with the bright siren analysis from <cit.> gives an H_0 measurement of 68.27^+4.20_-3.91 km/s/Mpc, also in agreement with the Legacy Survey result. Although the DELVE results present a slightly more precise measurement than those found with the Legacy Survey, the quality of the metrics of the photo-z measurements did not show the same robustness. The difference in reliability of photometric redshift measurements between DELVE and Legacy survey may be attributed to the absence of coadded data in DELVE, leading to lower signal-to-noise ratios (S/N) compared to the Legacy Survey catalog. Furthermore, Legacy offers uniform sky coverage, enabling analysis of events across both the Southern and Northern skies. Given that the precision of photometric redshifts primarily influences H_0 measurements, we selected the Legacy results as our primary outcome.


We also investigated the effect of adding galaxy fakes to the galaxy catalog for the S230919bj event. The presence of the galaxy fakes causes the peak of the posterior distribution to be slight shifted (a percent-level change) and the uncertainty is reduced by about 0.3%. 

Concerning the combined result of the 15 dark sirens, the presence of the galaxy fakes causes a 1% increase in precision, and the result combined with GW170817 remains unchanged, which shows the small impact that galaxy fakes have on the final conclusions of this work.

As other dark siren measurements, our H_0 measurement are valid under the Flat ΛCDM model and also presents a dependency with the background cosmology since the events go beyond z∼ 0.1 where changes to Ω_m and other cosmological parameters have a more significant impact on H_0 estimates. Therefore, we check if our results change with the choice of the Ω_m value within the  5σ interval found by the CMB measurement <cit.>. There is a minor shift, relative difference <4%, in the peak of the H_0 posterior distributions, and a relative difference <1% in the 68% CI. Therefore, the H_0 constraints presents here do not depend on the value of Ω_m, if it agrees with the Planck constraints.

The H_0 posteriors presented in this work were computed considering the full redshift PDF, computed through the deep learning algorithm described <ref>, for each of the possible host galaxies. The application of a photo-z PDF guarantees the use of a more reliable galaxies distribution first proposed in <cit.>, this treatment is different from other previous  <cit.> dark sirens works 
that assumed a Gaussian approach for these distributions. 
In Appendix <ref>, we discuss the comparison between the results of these different methodologies; the H_0 posteriors' behavior are slightly different, being more significant in closer events in which the effect of marginalization of thousands of galaxies is minimized.
As already noted by <cit.>, and also seen here, the use of the Gaussian approximation causes the distribution of galaxies dN/dz to be smoothed out, resulting in a flatter H_0 posterior. We note that the behavior of the H_0 posteriors generated by these two methodologies is slightly different, being more significant in closer events where the marginalization effect of thousands of galaxies is minimized. Despite these differences, the H_0 measurements are consistent with each other. 





































§ CONCLUSIONS


In this contribution, we use the data from the current best-localized and covered GW events from LVK O4a observing run to derive a dark siren measurement of the Hubble constant using the galaxy catalog method and precise photometric redshifts.  We obtained H_0=69.9^+13.3_-12.0 km s^-1 Mpc^-1, i.e. a ∼ 18% uncertainty on H_0 from dark sirens based on the catalog method alone and a total of 15 sirens. This is, at the best of our knowledge, an unprecedented precision for the catalog method and for the dark siren approach. 

We combine our result with 15 dark sirens with recent constraints over the one bright standard siren available GW170817, considering the constraints to the viewing angle from VLBI and the host galaxy peculiar velocity <cit.> and obtained H_0=68.00^+4.28_-3.85 km s^-1 Mpc^-1, a 6% measurement of H_0, reducing the previous constraint uncertainty in ∼ 10%. We note that the precision and uncertainty of the 15 dark sirens are similar to the GW170817 constraint without the viewing angle constraints from electromagnetic observations. It is worth noticing that current results are derived within the assumption of a Flat ΛCDM scenario.

Our current results emphasize that a combination of well-localized dark sirens and high-quality photometric redshifts can achieve a competitive H_0 constraint from gravitational waves, in particular, considering the absence of high confidence binary neutron star merger detections during O4a, the BNS merger rate can be in the lower-end of previous estimations. Therefore, the number of BBHs detections can be one order of magnitude higher than BNS. Furthermore, neutron star-black hole (NSBH) systems do not present substantial promise as multimessenger sources <cit.>. Nonetheless, it is prudent to acknowledge that these detection rates may significantly change in the near future, and the emergence of a singular, observable BNS event with an electromagnetic counterpart could potentially offer more compelling constraints than a dozen dark standard sirens. 




The current constraints from Dark Sirens by the catalog method only, achieve a precision of ∼ 18% and, in combination with Bright Sirens and additional constraints in the viewing angle, can achieve ∼ 6%. As the number of dark sirens events increase and we get closer to the level of statistical precision required to arbitrate the Hubble tension of ∼ 2% detailed studies to address potential systematics not included in this work should be carried out, especially considering different formation channels for BBH populations and catalog depth <cit.>.   




CRB acknowledges the financial support from CNPq (316072/2021-4) and from FAPERJ (grants 201.456/2022 and 210.330/2022) and the FINEP contract 01.22.0505.00 (ref. 1891/22). This material is based upon work supported by NSF Grant No. 2308193. AP thanks Constantina Nicolaou for sharing the GW170817 posterior. The authors made use of Sci-Mind servers machines developed by the CBPF AI LAB team and would like to thank A. Santos, P. Russano, and M. Portes de Albuquerque for all the support in infrastructure matters. 

The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; NSF's OIR Lab Proposal ID 2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NSF's OIR Lab Proposal ID 2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; NSF's OIR Lab Proposal ID 2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, The NSF's National Optical-Infrared Astronomy Research Laboratory (NSF's OIR Lab); the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NSF's OIR Lab. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham Nation.

The NSF's OIR Lab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science 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 University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, 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, Tecnologia e 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.

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).

The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.

The Legacy Surveys imaging of the DESI footprint is supported by the Director, Office of Science, Office of High Energy Physics of the U.S. Department of Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract; and by the U.S. National Science Foundation, Division of Astronomical Sciences under Contract No. AST-0950945 to NOAO.

The Photometric Redshifts for the Legacy Surveys (PRLS) catalog used in this paper was produced thanks to funding from the U.S. Department of Energy Office of Science, Office of High Energy Physics via grant DE-SC0007914.


 <cit.>,   <cit.>,   <cit.>

yahapj_twoauthor_arxiv_amp








§ GALAXY CATALOGS AND PHOTOMETRIC REDSHIFTS



In this work, we use catalog data from both DEcam Local Volume Exploration survey <cit.> and the Legacy Survey <cit.> as the main data used to select the galaxies within the LIGO/Virgo localization for each event. It is worth mentioning that the Legacy Survey catalogs incorporated publicly available DECam data, including DELVE footprint for the LS DR10.
The DELVE survey <cit.> observed a large fraction of the Southern sky covering an area of 21,000 deg^2, of which 17,000 deg^2 were homogeneously observed in the four broad bands (g, r, i, and  z) and photometric depth up to 24.3, 23.9, 23.5, and 22.8 for g,r, i, and z, respectively for 5σ point source detections. The DELVE catalogs comprise 2.5 billion sources, in which 618 million have data in all the four band. The Legacy survey explored a significant portion of the sky (∼ 33,500 sq deg) in grz bands, reaching depths of 24.0, 23.4, and 22.5 in grz for 5 σ detection.
The sky coverage of DELVE and Legacy surveys catalogs are shown in Fig. <ref>, together with the 50 and 90% CI of the GW events studied in previous works and the five new events explored in this study.


Photometric redshifts used in this work were determined using the same deep learning method presented in <cit.>. 
The fundamental component of the Deep Learning model comprises a neural network that analyzes tabular data, using a Legendre Memory Unit <cit.> with a Mixture Density Neural Network <cit.>.  Unlike traditional neural networks that provide single value estimations, MDNs outputs conditional probability densities through a linear combination of individual probability distributions (components), chosen to be Normal distributions in our case. This approach enables a more comprehensive characterization of predictions and errors assessment. 

The neural network output is a linear combination of C Gaussian kernels (𝒩(μ_i, σ_i), where {μ_i} is the mean and {σ_i} is the standard deviation) weighted by mixture coefficients {α_i}. We impose, for the mixture coefficients, that  ∑^C_i=1α_i=1 and 0<α_i<1. Therefore, the PDF can be written as


    PDF(z) = ∑^C_i=1α_i 𝒩(μ_i, σ_i).




This method was implemented to estimate PDFs for the DELVE DR2 and LS DR10 catalogs. The architecture used for the DELVE is detailed in <cit.>. The LS DR10[ We used the latest avaliable version of DR10, namely 10.1] where estimated using a similar architecture, with minor adjustments after a fine tunning considering the point statistics metrics such as reducing the width of the dense layers and reducing the number of components C in the mixture, from C=20 for DELVE to C=6 for Legacy. Input features for this process include griz magnitudes and color indices (g-r, g-i, g-z, r-i, r-z, and i-z). The spectroscopic data used as training set for DELVE catalogs is the same as the one presented in <cit.>, which was created from the crossmatch between DELVE catalog and several spectroscopic data available in different large sky surveys. For the Legacy survey we use the same training sample from the public photometric redshift catalog <cit.>.
 

The LS public data releases do not directly provides the apparent magnitude in photometric band. Therefore, we initially used the linear fluxes (columns ) to compute the magnitudes m by employing the conversion m=22.5-2.5log_10(f), and derived the magnitude errors also from the inverse variances of the fluxes (columns ) 
[https://www.legacysurvey.org/dr9/description/#photometryhttps://www.legacysurvey.org/dr9/description/photometry]. To mitigate star contamination, we adopted the same approach outlined by <cit.>, applying color cuts based on GAIA data, such as removing all known stars. Additionally, all the magnitudes 
were corrected for Milky Way extinction. Finally, we restrict our GW analysis to r band magnitudes to be lower than 21. 

We conducted a series of selection cuts based on the photometry quality and properties to ensure the utilization of galaxies with the best possible detections. We excluded all objects with unphysical colors, retaining only those satisfy the conditions:

    -1> g-r,  r-i, i-z <4







For the spectroscopic sample (refer to  for a detailed list of the spectroscopic catalogs also used in this work), we restricted our objects to 0.01<z_spec<1.5. After applying these cuts, our spectroscopic sample contains ∼ 2.2M galaxies.  We selected the training sample in order to have a uniform z_spec distribution of 0.01 to 1.0, and comprise all objects available with z_spec>1 (which are few in number), resulting in 550k (580k for DR9) galaxies for training the model. The same approach of uniform training was also used in Legacy Survey publicly available photo-z () to avoid possible systematic bias towards oversampled regions.
























We generate PDFs for galaxies in the test sample to validate the model by checking probabilistic and marginal calibration. Finally, we use the trained model to generate PDFs for galaxies in our target datasets. 
We explore the performance of the full PDF estimations by examining both the point-estimates photometric redshifts and the calibration of their PDFs. For a given galaxy, the value is defined as its respective PDF's most probable value (peak). Fig. <ref> presents the final photo-z distribution, found with the full photo-z PDF described in this appendix, in the 90% CI area of each gravitational wave event studied here. In order to highlight the overdensity regions, the photo-z distribution dN/dz was subtracted from the uniform distribution in comoving volume (dN/dz)_ com. 



To evaluate the quality of the point estimates, we employ the following metrics:


    
  * Median Bias: the bias, defined as Δ z = z_phot - z_spec, directly measures the deviation of our estimations from the target values (z_spec).

    
  * Scatter: the Normalized Median Absolute Deviation σ_NMAD, defined as
    
        
    σ_NMAD = 1.48×median(|Δ z - median(Δ z)/1+z_spec|),


    is a standard measurement of the bias scattering <cit.>. We aim for σ_NMAD to be as low as possible.
    The choice of σ_ NMAD instead of σ_68  is less sensitive to outliers. 

    
  * Outlier fraction: outliers are defined as objects which
    
        
    Δ z/1+z_spec > 0.15.

            We define η as being the fraction of outliers in any sub-sample of estimations. This definition of an outlier follows the same approach as adopted in <cit.> and <cit.>.
        


    All the aforementioned metrics were computed using the objects with spectroscopic correspondence inside the 90% confidence area of the events S231206cc, S230919bj, S230922g, GW191204_171526, GW190924_021846, GW200129_065458, GW200202_154313 and GW200311_115853. Additionally, with the exception of the odds constraints, the analysis refer only to the objects that satisfies the restrictions imposed to the events' catalogs. 

    
    

    Fig. <ref> presents the median bias, scatter, and outlier fraction, respectively, as a function of bins of 0.025 width. We computed the mean value of the metrics for each bin, generating one curve per event and then averaging these curves for all the aforementioned events. The solid lines represent the mean values across all the events, while the shaded regions represent the associated standard deviation. 

Considering the measurements for both DELVE and LS surveys, we obtain σ_NMAD and outlier fraction equal to 0.032 (0.016), 6.2(2.4)% for DELVE (Legacy), respectively. This results show that our measurements using Legacy data outperforms those obtained using DELVE data. We use the Legacy public photometric redshifts in Fig. <ref> catalogs for performance comparison in the point-like metrics. The LS public photometric redshifts were computed using the Random Forest (RF) method by <cit.>. In this method, they used r-band magnitude, g-r, r-z, z-W1, and W1-W2 colors as input features.
We compare our results with the public catalog from LS-DR10 as well as the public redshifts from LS-DR9 <cit.> utilized in <cit.>.

On the left panel of Fig. <ref> the MDN method exhibits an improvement of the averaged median bias in the 0.01<z<0.1 interval, compared both to the DR10 public available data and the data used in <cit.>, which is relevant since a considerable amount of the objects in the event skymaps lies in that interval. Conversely, for higher redshifts, both curves are compatible within the scatter. Throughout the entire range of , the bias for our method is very low and does not exceed 0.007. It is also consistent with bias 0 within the scatter until z∼ 0.4. However, it tends to be positive everywhere, indicating a weak cumulative overestimation effect. Additionally, in Fig.<ref> the scatter and outlier fraction align with publicly available redshifts values, thereby affirming the reliability of our method in generating accurate estimations. 
    
    
    
    
    
    

    Fig. <ref> illustrates our estimations plotted against the z_spec values available for all the events in this analysis. Despite the feature around  z=0.15, we observe well-distributed predictions around the diagonal line y=x. The averaged median bias for this entire sample is 0.002, while the averaged σ_NMAD and outlier fraction are 0.014 and 0.8%, indicating the robustness of our results. 
    
    



    
    
    
    
    
    
    

    
    
    
    
    
    

    
    


To validate the Probability Density Functions, we investigated the following metrics: 



    
  * PIT: 
    
    the Probability Integral Transform <cit.> represents the Cumulative Density Function of the PDF up to the z_spec value for each galaxy. We can assess the quality of the PDFs by examining the PIT distribution for a representative sample of galaxies. The distribution of PIT values is expected to be uniform between 0 and 1 <cit.>, indicating that the z_spec values can be considered as random events generated from independent PDFs. A slope in the PIT distribution indicates a bias in the estimations, while the concavity of the distribution reveals whether the PDFs are over- or under-dispersed <cit.>. 
    
    
  * Odds:
    
    the odds <cit.> measure the degree of confidence in the photo-z estimate derived from a given PDF. It is defined as the probability of the redshift lying within an interval around the photo-z value. For a given PDF, we can compute:
    
    
    Odds = ∫^z_phot+0.06_z_phot-0.06PDF(z)dz.

    
    The ideal distribution of odds for a galaxy sample well-represented by the training set should exhibit a pronounced peak near 1. However, it could also represent a distribution of under-dispersed PDFs. For this reason, we have to analyze the odds and PIT distribution simultaneously, in order to infer the quality of the PDFs.

    
  * Coverage diagnostic:  
    
    
    the coverage diagnostic (or High posterior density, HPD, diagnostic) stands as a well-established metric commonly employed to assess the quality of credible regions generated by simulation-based inference algorithms. Additionally, it provides a means to evaluate the accuracy of the estimated redshift distribution, as outlined in <cit.>. 
    The fundamental concept involves evaluating the probability that a specified credible region within the inferred distribution contains the true value. This assessment provides insights into whether the estimated distribution is overconfident, calibrated, or underconfident <cit.>. The coverage diagnostic was executed by selecting the pair (spec-z, PDF) and sampling 1000 values from the photo-z's PDF, thereby generating a frequency distribution with 1000 bins spanning the range from 0 to 1. 
    
    
    
    First, a credible region is defined for the estimated distribution using the Highest Density Regions (HDR). This region represents the smallest area that contains at least 100(1 - α)% of the mass of the inferred photo-z distribution, establishing an interval for a given credibility level (1-α ). The expected coverage is the frequency with which the true parameter (spec-z) value falls within this highest density region, essentially indicating how often it falls inside the calculated interval.  
    If our model produces well calibrated distributions we expect the spec-z value should be contained inside the calculated interval from the (1-α) HPD region of the estimated distribution exactly (1 - α)% of the time.
    If the coverage probability is less than the (1 - α) credibility level is the sign of under-estimation of the PDF's variance, and it could lead to unreliable approximations since it excludes physical values of . Conversely, if the coverage probability is larger than the (1 - α) credibility level, then this indicates that the estimated PDFs are over-estimating their variance, in average.
    

    



    


 The middle panel of Fig. <ref> show the odds distribution for the same set of galaxies present on the point-estimate analysis, while the left panel illustrates the PIT distribution for a subset of these objects satisfying the constraint odds>0.7. 
The latter exhibits a marginally concavity at lower values of PIT and a negative slope throughout higher values. Although not pronounced, this features collectively suggest a tendency towards positive bias, implying a slight overestimation of the most probable 's values. 
    
    As discussed above, there is trending in the median bias curve (Fig. <ref>) of overestimating . This overestimation effect is further illustrated in Figure <ref>, particularly evident around equal to 0.48. 
    
    Despite this anomaly at the higher redshifts, the scattering of remains well behaved, as expected. Consequently, the inclination of the PIT distribution might be attributed to a localized systematic bias rather than indicating any inherent global bias within the model. Although this could impact the analysis, it is worth stressing that our Bayesian analysis performs a bias correction based on this curve in the same manner as previous work by <cit.>. A similar shape of PIT distribution was also encountered by <cit.> using MDN to estimate redshifts from images in the Sloan Digital Sky Survey Data Release 9.
    
    
    
    The odds distribution shows a gradual increase in frequency with the odds values. This result suggest the absence of a systematically overconfident or overdispersed PDFs, indicating a high level of confidence in the estimated PDFs. From Fig.<ref>, a perfectly calibrated estimated distribution aligns with a diagonal line (dashed black line), indicating an expected coverage probability matching the credibility level.  We notice that the expected coverage curve matches the diagonal, indicating that our model produces well calibrated PDFs. This linear relation confirms the reliability of our methodology. 

    In summary, our methodology estimates with an accuracy comparable to those publicly available in LS-DR10, while also generating well-calibrated PDFs. It is noteworthy that our approach can predict redshifts in objects with detections in all g,r,i,z bands, and in any combination of three bands containing detections in g band.
    
   


















































§ METHOD

Here, we outline the statistical methodology known as the dark siren approach, initially introduced in <cit.> and subsequently modified by <cit.>.
The method consists of using the Bayesian formalism to infer the Hubble constant (H_0) parameter through the gravitational wave detection data d_ GW and the d_ EM set of photo-z measurements of the possible host galaxies made using the deep learning algorithm (see Appendix <ref>). As the GW and EM measurements are done independently, the joint GW and EM likelihood can be defined as the product of the two individual likelihood, p(d_ GW, d_ EM| H_0) = p(d_ GW|H_0) p(d_ EM|H_0). 

From the Bayesian framework, the H_0 posterior of one GW event can be written in the final form as


    p(H_0|d_ GW, d_ EM) ∝p(H_0)/β (H_0)∑_i 1/𝒵_i∫ z  Δ z   p(d_ GW|d_L(z,H_0),Ω̂_i) p_i(d_ EM|z,Δ z)  p(Δ z) r^2(z)/H(z)  .


where p(H_0) is prior on H_0 which we assume to be flat over the range [20,140] km s^-1 Mpc^-1, β (H_0) is the normalization factor that describes the selection effects in the measurement process, r(z) is the comoving distance, H(z)=H_0(Ω_m(1+z)^3+1-Ω_m)^1/2 is the Hubble parameter in a Flat ΛCDM model, 𝒵_i = ∫ p(d_ EM | z_i) r^2(z_i)/H(z_i)   z_i are evidence terms that normalize the posterior, and p(Δ z) is the prior on the photo-z bias that we measure from our photo-z validation sample (see Fig. <ref> in appendix <ref>). In the following paragraphs we will detail how the function β(H_0) was computed. The GW and EM likelihoods are written as, respectively:

    p(d_ GW|d_L(z,H_0), Ω̂_i)∝p(Ω̂_i)N(Ω̂_i)/√(2π)σ(Ω̂_i)exp[-(d_L-μ(Ω̂_i))^2/2σ^2(Ω̂_i)], 
    
    p(d_ EM|z, Δ z) = ∏_i p(z_i|z, Δ z),

where we explicitly consider the dependence of cosmological parameters on the luminosity distance, as measured by GW, and the solid angle Ω̂_i corresponding to each observed galaxy i. Following <cit.>, the GW likelihood is approximated by a Gaussian function. The EM likelihood is the product of the probability distribution function of the photometric redshift z_k for each k galaxy, where we consider the correction of the shifted in redshift for the photo-z biases, Δ z, in the data. This phenomenon arises due to the absence of a uniform distribution in redshift or color beyond a certain magnitude or color selection threshold, which makes the deep learning algorithm begin to oversample the redshift around the distribution peaks, thereby introducing systematic biases. The individual H_0 posterior distributions found in this work are presented (in colors) in Fig. <ref>.

 

The selection effects were computed based on the criteria outlined in <cit.> and <cit.>. This involves the joint gravitational wave-electromagnetic likelihood, which is marginalized over all conceivable GW and EM data. Assuming that the events are isotropically distributed on a large scale, this term can be written in its compact form as:


    β = ∫_0^z_ max p^ GW_ sel(d_L(z,H_0))p(z)dz


where p_ sel^ GW is the probability of a GW event at a given luminosity distance d_L be detected, p(z) is the galaxy catalog distribution, and z_ max is the maximum true redshift at which we can detect the host galaxies. The function p_ sel^ GW is equal to one if all events in a given redshift z satisfy the detection condition, i.e. the detector network SNR is >12 and the localization volume satisfies our selection criteria (A_50%<1000 ^2 and d_L<1500 Mpc), and zero if none of the events located in a redshift z satisfies those conditions. Therefore, p_ sel^ GW function is, in sum, an efficiency curve, i.e. a smooth function that falls from 1 to 0 over the redshift range for each H_0 value.

The selection function was modeled by simulating 30,000 BBH mergers using  software <cit.> with the frequency domain approximant IMRPhenomD. We simulate the BBH mergers assuming that they follow a power law plus peak mass distribution with the same parameters as described in <cit.>, the spins distributions follow a uniform distribution between (-1,1), and a uniform distribution in comoving volume assuming a dependency with the merger rate evolution described by the Madau-Dickinson star formation rate <cit.> and fixing the Planck 2018 cosmology for twenty different H_0 values within our prior range. We assumed that the Legacy galaxies catalog is complete up to the known absolute magnitudes for each of the GW events in the analysis, which correspond to a maximum redshift of z_ max=0.65. All the 30,000 injected events passed by the matched-filter analysis, where we computed the signal-noise-ratio (SNR) assuming the O4 sensibility curve as published by LIGO in Document P1200087 (<https://dcc.ligo.org/LIGO-T2000012/public>) for the Hanford and Livingston LIGO network detectors. The network SNR defines the detection condition above 12 and at least 2 detectors have a single–detector SNR above 4. A Gaussian noise was added in all the measurements. Lastly, the  skymaps were reconstructed for each detection event, where we assume a luminosity distance that follows ∝ d_L^2.

The last selection cut that we should consider is that events serving as dark sirens have at least 70% of their 90% CI comoving volume covered by the Legacy survey. As noted by <cit.>, this selection cut can be ignored since the GW antenna pattern is not correlated to the survey sky footprint and we do not expect an H_0 dependence. Here, we also ignore this selection effect. 




§ COMPARISON WITH H_0 CONSTRAINTS FROM POINT ESTIMATES OF PHOTOMETRIC REDSHIFTS


During this work, the H_0 posterior calculation was performed using the full photo-z PDFs estimated through the deep learning method, described in Appendix <ref>. In this appendix we will analyze the effect on the H_0 posterior originated by the choice of the galaxies photo-z PDFs, the  photo-z PDF estimated by the MDN technique described in Appendix <ref>, which we will call the full photo-z PDF, or a Gaussian approximation, where we will approximate each full photo-z PDFs to a Gaussian function whose mean coincide with the peak value and the value of the standard deviation is equal to that found with the full photo-z PDF.

In the sense that the Gaussian approximations were constructed with the same peak and standard deviation as that found with the PDFs estimated with the MDN. We include the results of the public Legacy [The photo-z were computed using the Random Forest technique with data from the Legacy DR10.1, these results are available in <https://www.legacysurvey.org/dr10>.] in the figs. The top panel of Fig. <ref> shows the posterior of H_0 found considering the full photo-z PDF (blue solid curve) and the Gaussian approximation (black dashed curve for the full PDF, this work, and the Legacy DR10.1 public results in gray dashed curve), and the residuals curves are shown on the bottom panels. Comparing the results of the posterior of H_0 generated considering the Gaussian approach (black dashed curve) and the full photo-z PDF, we note that the results are in agreement at the percentage level for the more distant events, where the effect is suppressed by marginalization over a larger number of galaxies due to its larger localization volume. The effect of the use of the Gaussian approximation is more evident for the events S230627c, GW200129 and GW191204, which present a discrepancy of ≈ 4%, 16%, and 40%, respectively, in the peak distribution and can reach values >13% at the ends. As expected, the Gaussian approximation makes the peak of the H_0 posterior wider, implying a less precise H_0 measurement. This result is a consequence of the smoothing of overdensity regions in the photo-z distribution. Although the results do not indicate that the choice of photo-z PDF leads to disagreement H_0 measurements, future dark siren measurements can reach a level of precision where the differences between these results may be statistically significant. We can also see that the Legacy public results leads to, in most events, a less restrictive H_0 posterior compared to those achieved through the full photo-z PDFs, leading to a increase of up to 6% in the H_0 uncertainty.

As shown by <cit.>, the choice for a wider redshift cut causes galaxies to be added in deeper redshift, whose photo-z PDF has the a significant Gaussian tails at high redshift, which implies a increase of the H_0 posterior at high H_0. However, the full photo-z PDF is able to reduce (the residual values, between the posterior considered a LIGO/Virgo luminosity distance posterior of 90% and 99.7%, at the high-H_0 end reduced by ≈ 25% compared to the result achieved with the Gaussian approximation) this dependence of the H_0 posterior behavior with redshift cut, which evidences its advantage against the Gaussian approximation in performing a measure of H_0 free of any systematic imposed by assumptions made in the methodology.  



In Fig.<ref> we present the photo-z PDFs estimated by the deep Learning method (dark blue curve) and the Gaussian approximation (dashed curve in light blue) for galaxies at different redshifts. We notice that the full photo-z PDF has a tendency to be narrower than the Gaussian approximation, which favors a more precise measurement of H_0. In addition, we see that the difference between the estimated PDF and the Gaussian approximation becomes more evident in higher redshifts. This highlights the importance of the MDN technique mainly for galaxies located in deep redshifts, that are already used for the LIGO/Virgo GW dark sirens and will be increasingly indispensable for the coming years where it is expected to observe increasingly distant GW events. We can also observe that the full photo-z PDFs become increasingly non-Gaussian in deeper redshifts. In order to quantify this non-Gaussianity in Fig. <ref> we present the relation between the kurtosis (fourth standardized moment μ_4, that is presented in the left panel) and the skewness (third standardized moment μ_3 shown in the right panel) of the full PDFs with z_ phot, where we compute the mean (solid line) and the standard deviation (shaded region) of this quantities in photo-z bins of size 0.05. The values of kurtosis and skewness are distinct from the values of a Gaussian distribution (kurtosis and skewness equals to 3 and 0, respectively) in almost the entire z_ phot interval. We note a tendency of kurtosis values to be lower with increasing z_ phot, which indicates the enlargement of PDFs implying in less precise z_ phot measures. This behavior describes the already expected difficulty of precise measurements of distant galaxies. The skewness results indicate that the full PDFs are right-skewed in low z_ phot, being reduced with the increase of z_ phot until reaching the value of zero (indicating a symmetric PDF) in z_ phot≈0.5, in which from this value begins to be left-skewed. 

 

 




