The detection of gravitational waves GWs with an electromagnetic counterpart enabled the first Hubble constant inlineformulammlmath displayinlinemmlmrowmmlmsubmmlmrowmmlmiHmmlmimmlmrowmmlmrowmmlmn0mmlmnmmlmrowmmlmsubmmlmrowmmlmathinlineformula measurement through the standard siren method. Current constraints suggest that inlineformulammlmath displayinlinemmlmommlmommlmn20mmlmnmmlmimmlmimmlmn30mmlmnmmlmommlmommlmathinlineformula of LIGOVirgoKAGRA LVK binary black hole BBH mergers might occur in active galactic nuclei AGN disks. The claim for a possible association of several BBH mergers with flaring AGNs suggests that cosmological analyses using BBHs and AGNs might be promising. We explore standard siren analyses through a method that takes into account the presence of background flaring AGNs without requiring a unique host galaxy identification and apply it to realistic GW simulations. Depending on the fraction of LVK BBHs that induce flares we expect to constrain inlineformulammlmath displayinlinemmlmsubmmlmiHmmlmimmlmn0mmlmnmmlmsubmmlmathinlineformula at the inlineformulammlmath displayinlinemmlmommlmommlmn6mmlmnmmlmimmlmimmlmn7mmlmnmmlmommlmommlmathinlineformula inlineformulammlmath displayinlinemmlmommlmommlmn4mmlmnmmlmimmlmimmlmn5mmlmnmmlmommlmommlmathinlineformula precision with inlineformulammlmath displayinlinemmlmommlmommlmn2mmlmnmmlmimmlmimmlmn3mmlmnmmlmtext mmlmtextmmlmtext mmlmtextmmlmiyearsmmlmimmlmathinlineformula or inlineformulammlmath displayinlinemmlmommlmommlmn160mmlmnmmlmimmlmimmlmn240mmlmnmmlmathinlineformula events inlineformulammlmath displayinlinemmlmommlmommlmn1mmlmnmmlmtext mmlmtextmmlmtext mmlmtextmmlmiyearmmlmimmlmathinlineformula or 500 events of LVK at design inlineformulammlmath displayinlinemmlmrowmmlmi mathvariantnormalAmmlmimmlmommlmommlmrowmmlmathinlineformula sensitivity assuming that systematic BBH followup searches are performed. We also show that in a scenario where only inlineformulammlmath displayinlinemmlmommlmommlmn1mmlmnmmlmommlmommlmathinlineformula of the BBHs induce detectable flares it is possible to achieve an inlineformulammlmath displayinlinemmlmsubmmlmiHmmlmimmlmn0mmlmnmmlmsubmmlmathinlineformula precision from 7.5 to 15 with inlineformulammlmath displayinlinemmlmommlmommlmn240mmlmnmmlmathinlineformula events. In addition a inlineformulammlmath displayinlinemmlmommlmommlmn5mmlmnmmlmimmlmimmlmn27mmlmnmmlmommlmommlmathinlineformula precision is achievable with complete AGN catalogs and 1 year of LVK run without the need of any flare identification.
['Astrophysics and astroparticle physics', 'Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - High Energy Astrophysical Phenomena']
Standard siren cosmology with gravitational waves from binary black hole mergers in active galactic nuclei
{'Standard Siren Cosmology with Gravitational Waves from Binary Black Hole Mergers in Active Galaxy Nuclei': 'C. R. Bom 1, 2, ∗ and A. Palmese 3, †  \n1 Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, 22290-180 Rio de Janeiro, RJ, Brazil 2 Centro Federal de Educação Tecnológica Celso Suckow da Fonseca, Rodovia Márcio Covas, lote J2, quadra J - Itaguaí (Brazil) 3 McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA  \nThe detection of gravitational waves (GW) with an electromagnetic counterpart enabled the first Hubble Constant H 0 measurement through the standard siren method. Current constraints suggest that ∼ 20 -30% of LIGO/Virgo/KAGRA (LVK) Binary Black Hole (BBH) mergers might occur in Active Galactic Nuclei (AGN) disks. The claim for a possible association of several BBH mergers with flaring AGNs suggests that cosmological analyses using BBH and AGNs might be promising. We explore standard siren analyses through a method that takes into account the presence of background flaring AGNs, without requiring a unique host galaxy identification, and apply it to realistic GW simulations. Depending on the fraction of LVK BBHs that induce flares, we expect to constrain H 0 at the ∼ 6 -7% ( ∼ 4 -5% ) precision with ∼ 2 -3 years or ∼ 160 -240 events ( ∼ 1 year or 500 events) of LVK at design (A+) sensitivity, assuming that systematic BBH follow-up searches are performed. We also show that in a scenario where only ∼ 1% of the BBHs induce detectable flares it is possible to achieve an H 0 precision from 7 . 5% to 15% with ∼ 240 events. In addition, a ∼ 5 -27% precision is achievable with complete AGN catalogs and 1 year of LVK run, without the need of any flare identification.', 'I. INTRODUCTION': "The detection of gravitational Waves (GW) by LIGO and Virgo has enabled a variety of novel measurements of cosmological parameters. The 'standard siren' method, first proposed by [1], relies on a luminosity distance measurement from the GW detection, and together with an independent measurement of the redshift of the GW event, it allows us to probe the expansion of the Universe through the distance-redshift relation. The current 4 -6 σ tension between the measurement of the Hubble constant H 0 from Cosmic Microwave Background [2] and Cepheid-anchored Type Ia Supernova [3] analyses makes an independent, standard siren measurement of H 0 from gravitational waves a promising endeavor to shed light on the origin of this tension through an independent probe. However, despite enormous efforts by the astronomical community, only one optical counterpart to a GW event, GW170817 [4, 5], has been confidently identified, and the identification of its host galaxy along with its redshift enabled standard siren measurements of H 0 [6-8]. About ∼ 50 similar associations are needed to produce a measurement at the precision level needed to discern between the two H 0 measurements currently in major discrepancy [9]. In the cases where no counterpart is identified it is possible to use a different approach, the 'dark siren' method, by taking advantage of assumptions on the mass distribution of compact objects (e.g. [10, 11]), or by taking into account the ensemble of potential host galaxies  \n[1, 12, 13]. Recent dark siren analyses [14-18] reach a precision on H 0 down to ∼ 20% . Hundreds to thousands of such events (or a few 'special' events [19]) will be required to reach a 2% uncertainty on H 0 (e.g. [20]). For current generation GW detectors, most of the dark siren works make use of binary black hole (BBH) mergers. BBHs are usually expected to originate as isolated binaries (e.g. [21]) or through dynamical channels (e.g. [22-27]). A promising site for the formation of massive stellar mass binaries through hierarchical mergers is in Active Galactic Nuclei (AGN) [23]. Recent works [28-30] show that ∼ 20% of the BBHs detected by LIGO/Virgo are expected to occur in AGNs. In the AGN environment, the electromagnetic signature associated with BBHs is contingent upon the interplay between the merger and the adjacent disk gas. Generally, a BBH within an AGN disk is displaced from the point of merger, exhibiting a recoil kick velocity determined by the mass asymmetry and spin orientation of the progenitor system. The gas in the disk around the progenitor BBH tends to follow the merged remnant, yet it encounters the neighboring AGN disk gas, consequently generating a luminous shockwave within the disk. Gas accretion onto the kicked black hole may give rise to a Bondi tail, or a jet may be launched and produce an afterglow. These phenomena may be observable in a wide range of electromagnetic radiation wavelengths, including in the optical/ultraviolet [31-35].  \nSeveral AGN flares have been flagged as possible electromagnetic counterparts [36, 37], but given the variable nature of AGNs, it is challenging to establish a confident association with low probability of chance coincidence for GW events with a broad sky localization, and multiple flares may occur within the GW region around the same time. On the other hand, one could consider using all of  \nthe AGNs that exist within the GW localization, in case no follow-up observations exist or a flare is not expected (e.g. for obscured AGNs). Even though it is challenging to identify which flare or which AGN is indeed associated with a single GW event, after enough observations of GW events along with AGN flares or AGN catalogs, it will be possible, statistically, to say whether the GWAGN associations are consistent with a background of chance coincidence observations or they are indeed related, given that the spatial and temporal distributions of chance associations will be different from that of a real association [38-40]. If the association is deemed significant, then the AGN population also provides the redshift information needed for standard siren measurements.  \nWe explore whether making use of the population of BBH occurring in AGNs can improve standard siren constraints, considering that AGNs are a small subsample of the generic galaxy population such as that used in dark siren analyses so far, and that the number of host galaxies one marginalizes over dictates the final cosmological parameter precision. Furthermore, the distance reach of BBH mergers (luminosity distance ∼ 1400 -1600 Mpc for the LIGO detectors in the current observing run [41]), as opposed to that of binary neutron star (BNS) mergers ( ∼ 160 -190 Mpc), will allow us to eventually probe cosmological parameters beyond H 0 , such as the Universe matter density and the dark energy equation of state. Note that this approach is different from that of previous work assuming a unique host galaxy is identified with high confidence [42-45] through identification of a unique, confirmed electromagnetic counterpart. Motivated by the findings of [37] and the expected sky localizations of GW events in the coming years, it is reasonable to assume that at least for current generation GW detectors high confidence associations on single events will be challenging. We base our forecast on the assumption that the associations can only be made at a population level, and that is also the way our standard siren approach is implemented. Because of this, our constraints will be very different from those from previous works. First, precision: because the uncertainty in counterpart association is taken into account here, we expect our forecasts to be less stringent and more realistic than a forecast assuming a high confidence association. Second, accuracy: our simulations include background flares as potential counterparts, yet our method yields unbiased constraints, as we show in the following sections. On the other hand, if one uses the classical bright standard siren approach with a flare that is unlikely to be associated, it is clear that the analysis may yield biased constraints, since one would use an incorrect value of the redshift.  \nIn this work, we explore the promising avenue of taking advantage of the BBH events that merge in AGN disks to constrain cosmological parameters. The formalism was first presented in [39], where the authors show how to constrain the fraction of LIGO/Virgo BBH that give rise to AGN flares, while also deriving cosmological parameters. The method is generic and can be used to constrain  \nthe fraction of BBHs hosted by AGNs (as also explored in [38]) and cosmological parameters without the requirement of a flare, hence of follow-up observations, assuming that a complete AGN catalog is available.", 'A. Standard Siren Cosmology with Noisy Source Identification Environment': "The association between BBHs and AGN flares might be addressed as a signal-to-background distinction, as done for the formalism presented in [39]. The total number of AGN flares per solid angle Ω and redshift z , dN d Ω dz , can be construed as a composite model consisting of an AGN flare associated to the gravitational wave (GW) event, if any, and the background number density of AGN flares within a specified temporal interval. dN d Ω dz is dependent on λ , the fraction of BBHs inducing AGN flares. For a set of events, we can then write down a posterior to constrain λ .  \nThe galaxy catalog approach of the dark siren method uses the redshift information of possible host galaxies of the GW event and the luminosity distance obtained from GW detectors. Differently from a traditional dark siren approach, we consider AGN flares within the high probability sky location of the GW event occurring within a time window from the BBH event detection, to constrain the fraction of events that induce a 'signal' flare, i.e. λ . Note that by including this factor and letting it be λ < 1 , we are allowing for the fact that not all BBH, but only a fraction, may occur in AGN. Similarly, one can consider known AGNs (instead of flaring AGNs or all possible host galaxies) where a follow-up is not available, in which case we refer to λ AGN . Thus, the proposed method takes advantage of an intermediate perspective between the dark and bright siren approaches considered so far.  \nLet N be the number of GW events with GW data { x GW i } N i =1 ≡ { Ω GW i , d GW L,i } N i =1 and with AGN data { x AGN i } N i =1 ≡ { { Ω AGN ij , z AGN ij } k j =1 } N i =1 , where Ω is a sky line of sight, d L is the luminosity distance and z is the redshift. The subscript i refers to a GW event, while j refers to an AGN potentially associated to the event i . After applying Bayes' theorem and marginalizing over the fraction of GW events hosted by AGNs λ , the posterior on the cosmological parameters ⃗ θ is given by:  \np ( ⃗ θ | { x AGN i } N i =1 ) , { x GW i } N i =1 ) ∝ p ( ⃗ θ ) ∫ dλp ( λ ) N ∏ i L i ( ⃗ θ, λ ) (1)  \nThe likelihood of observing the GW data x GW i and k AGNs or AGN flares with sky positions and redshifts { Ω AGN ij , z AGN ij } k j =1 , is given by an inhomogeneous Poisson process. We then marginalize over the position of  \nthe GW source (Ω GW i , z GW i ) , to obtain the likelihood L i ( ⃗ θ, λ ) as:  \n(2)  \nL i ( ⃗ θ, λ ) ∝ ∝ k ∏ j =1 [ λp ( x GW i | Ω AGN ij , d L ( z AGN ij , ⃗ θ )) p 0 (Ω AGN ij , z AGN ij ) + R B (Ω AGN ij , z AGN ij , ⃗ θ ) ] e -µ i .  \nIn Eq. 2 the term that multiplies λ is effectively the GW posterior on the sky location and redshift evaluated at the AGN position and redshift, for a given cosmology. The second term in brackets is the background rate, and it is defined as R B = T dB d Ω dzdt , i.e. the rate of background flares or AGNs per unit solid angle and redshift within a time period T . At last, the term µ i is the expected number of observed AGN flares or AGNs:  \nµ i ≡ ∫ dN i d Ω dz P AGN det (Ω , z ) d Ω dz , (3)  \nwhere P AGN det is the detection probability of the AGNs/flares. For a full derivation of this likelihood see [39]. We sample the posterior on the cosmological parameters and λ using the python package emcee .  \nWe use realistic simulations of GW events for LIGO/Virgo/KAGRA (LVK) close to the expected design sensitivity (and fourth observing run, O4, sensitivity) and for their A+ sensitivity, along with their expected 3D localization reconstructions, produced as described in Section II B. To estimate the number of AGNs within a GW localization volume, we use a fiducial number of n AGN = 10 -4 . 75 Mpc -3 [46-48], and also consider a second scenario with higher background contamination ( n AGN = 10 -4 Mpc -3 , motivated by the fact that BBH merger may extend to low accretion rates AGNs [24]). The number of background AGNs flaring depends on several factors. As most of AGNs present some variability, the definition might depend on the rise in brightness timescale τ , and also be sensitive to the variability model chosen. For instance, the first flare claimed to be of BBH origin in AGN J124942.3+344929 presented a a rise of ∆ m ∼ 0 . 4 mag over ∼ 50 days [36]. For this event, [39] used structure functions (SF) to derive a probability of flare α flare ∼ 10 -4 with ∆ m > 0 . 4 . On the other hand [36] estimated, α flare = 5 × 10 -6 considering a fit in the sample of AGNs data from the Zwicky Transient Facility (ZTF). As the number of background flares can be influenced by the exact theoretical modeling, which is yet to be well understood for BBH events, and the timescales probed by a specific follow-up campaign, we choose the more conservative estimation from [39] (i.e. higher background flare contamination than that estimated in e.g. [36]) scenario, one similar to the AGN J124942.3+344929 flare background estimate, α flare = 10 -4 .  \nEven though [28] find that up to ∼ 80% of BBHs detected by LIGO/Virgo/KAGRA might originate from the AGN channel, many could be hosted in an accretion disk  \nobscured from view (type-II AGNs). Although flares occurring in this type of AGN might still become detectable at later times and at longer, IR wavelengths (see, for instance [49]), this may impact the fraction of BBHs inducing flares, hence the λ we consider in the case of flaring AGNs should be lower than that of [28], which will include some obscured AGNs. Using the prescription from [50] the fraction of unobscured AGNs varies with the bolometric luminosity, it can be ∼ 20% ( ∼ 26% ) for AGNs of luminosity 10 44 ( 10 46 ) erg s -1 visible at optical wavelengths. In the case of GW190521 in association with AGN J124942.3 + 344929, considering AGNs with g < 20 . 5 , the type-I AGNs comprise ∼ 22% of the entire AGN sample [39]. Thus, if only BBHs flares from type-I AGNs were detectable, this would limit λ from < 0 . 8 to < 0 . 18 . Therefore, in our fiducial scenario, we restrict the analysis to λ ∈ [0 . 05 , 0 . 3] , where the lower limit is chosen by selecting type-I AGNs for ∼ 20% of BBH in AGNs found by [30]. Furthermore, since [29] find a lower limit for the AGN channel of ∼ 3% of all BBH events, we also consider an alternative, more pessimistic scenario, where λ can go down to 1/4 of the lower bound in our fiducial scenario ( λ ∈ [0 . 0125 , 0 . 025] ). We note that the lower bound here is not ∼ 22% of the 3% value from [29], as bringing λ below 1% would not produce any interesting results, at least for the upcoming gravitational wave runs.", 'B. Gravitational Wave simulations': "We analyze a dataset of simulated BBH mergers generated according to the expected observations from the current LIGO-Virgo-KAGRA (LVK) observing run (the fourth run, 'O4'), and the following LVK A+ upgrade (potentially close to what expected for the fifth observing run, O5). Our methodology is similar to that employed in [39, 51, 52]. Gravitational wave events are generated using the BAYESTAR [53, 54] software, which makes use of tools from LALSuite [55]. We adopt sensitivity curves expected for Advanced LIGO and Virgo at design sensitivity and for KAGRA at a BNS inspiral range of 80 Mpc, and at A+ sensitivity, all taken from https: //dcc.ligo.org/LIGO-T2000012-v1/public , considering a duty cycle of 70% [41]. In the A+ scenario, we also include a LIGO-India detector. The BBH population follows a mass distribution in the form of a power law plus peak with parameters from [56]. We set the input cosmological parameters based on a flat Λ CDM cosmology with parameters from [2]. In Figure 1 we show the 90% CI area versus luminosity distance (top panel), and the comoving volume versus fractional distance error (bottom panel) for the simulations used in this work. Clearly, with A+ we will be able to detect more events and up to larger distances than with LVK at design sensitivity, resulting in a larger fraction of high comoving volume events. At the same time, with A+ a large number of extremely well localized 10 sq. deg. events will be  \nFIG. 1. The LVK (O4, indigo) and LVK A+ (O5, gray) parameters set obtained from simulations of binary black hole mergers. Top: distance and 90% credible interval sky localization area of our simulated events. Bottom: The GW 90% credible interval Volume and the fractional error in the GW Luminosity distance.  \n<!-- image -->  \navailable.  \nWe use the events' 90% comoving volume to compute the expectation value of the number of background flares or AGNs for each event, and draw for each simulated event a random number of background objects following a Poisson distribution. The background objects follow a uniform in comoving volume redshift distribution, and an isotropic distribution in the sky area (i.e. we ignore AGN clustering for the purpose of this work). On the other hand, whether a GW event has an associated 'signal' flare, depends on the input value of λ , again following a Poisson distribution. Signal flares follow a different sky and redshift distribution, which is dictated by the GW event posteriors.", 'A. AGN flares case': "<!-- image -->  \nFIG. 2. Hubble constant precision as a function of number of events from LVK (similarly to an O4 design sensitivity) and LVK A+ (similar to O5) for different input values of the fraction of BBHs inducing AGN flares λ . The upper dot-dashed horizontal line presents the expected H 0 precision from 250 dark sirens using BNSs from [57], considering a 50% complete galaxy catalog and unweighted galaxies. The lower dot-dashed line presents the forecast from [58] after ∼ 50 BNS bright sirens. The shaded region shows the uncertainties over 100 realizations expected for a golden dark siren where we can identify the host from HLV+ with spectroscopic redshift [19]. The top (bottom) panel considers a flat prior over Ω m ∈ [0 , 1 . 0] ( [0 . 25 , 0 . 35] ). The number density of AGN is n AGN = 10 -4 . 75 Mpc -3 for the filled symbols and 10 -4 Mpc -3 for the empty symbols.  \n<!-- image -->  \nIn Figure 2, we present the precision on H 0 obtained in a flat Λ CDM cosmology with two choices of uniform priors for the matter density Ω m , Ω m ∈ [0 . 0 , 1 . 0] in the top panel and Ω m ∈ [0 . 25 , 0 . 35] , corresponding to ∼ 5 σ  \n<!-- image -->  \nFIG. 3. Hubble constant precision as a function of λ from LVK design and LVK A+ for the most pessimistic scenarios considered here, with a lower limit on the fraction of BBH producing flaring AGNs of ∼ 10 -2 . The upper dot-dashed horizontal line presents the expected H 0 precision from 250 Dark Siren using BNSs from [57], considering a 50% complete galaxy catalog and unweighted galaxies. The lower dotdashed line presents the forecast from [58] after ∼ 50 BNS bright sirens. The shaded region shows the uncertainties over 100 realizations expected for a golden dark siren where we can identify the host from HLV+ with spectroscopic redshift. The number density AGN is n AGN = 10 -4 . 75 Mpc -3 for the filled symbols and 10 -4 Mpc -3 for the empty symbols. These results consider a flat prior on Ω m ∈ [0 . 25 , 0 . 35] .  \n<!-- image -->  \naround the best fit Planck Cosmology value from [2], in the bottom panel. The prior in the Hubble constant is uniform H 0 ∈ [20 , 140] km/s/Mpc, and λ ∈ [0 . 0 , 1 . 0] , as a function of the number of events in the LVK design and LVK A+ scenarios, for different input values of λ , and for the two different AGN number density considered ( n AGN = 10 -4 . 75 Mpc -3 for the filled symbols, and n AGN = 10 -4 Mpc -3 for the empty symbols of Figure 2). The points shown represent the median and the errorbars represent the 68% credible interval using 10 different runs. While the number of runs considered is on the low side for a sample that should ideally be used to derive statistics, the computational time required for a larger number of runs for all the scenarios considered would be prohibitive.  \nThe number of events shown corresponds to the median expectation in [41] per year in O4, i.e. ∼ 80 events for LVK and the median expectation in [51] of ∼ 500 events in LVK A+ (since O5 prospects are not included in the official LVK forecasts of [41]). O4 is currently scheduled for ∼ 22 months of data acquisition, hence the median in the whole O4 run could be pushed close to ∼ 150 . Therefore, for the current observing run, the expectations might be close to what we report for 160 events. Note that [51] find a median of 106 BBH detections in O4, and considering the uncertainties we could  \nFIG. 4. Posteriors from 1 year (or 500 events) of LVK A+ within a flaring AGN scenario with n AGN = 10 -4 . 75 Mpc -3 and α flare = 10 -4 . Top: H 0 posterior considering a Λ CDM model for different input values of λ . The solid lines represent the result for a restricted flat prior of Ω m ∈ [0 . 25 , 0 . 35] with the precision quoted in the legend. The dashed lines assume a flat wide prior Ω m ∈ [0 . 0 , 1 . 0] , and their respective recovered precision is quoted in parentheses in the legend. The vertical line is the input value of H 0 . Bottom: The dark energy equation of state parameter w posterior considering a w CDM cosmology for different input values of λ . The legend presents the median w and the best fit (in parentheses), with the 68% credible intervals. The dot-dashed lines present the 68% credible intervals and the mode. The input value of w in our simulations is -1. The gray-shaded region shows the forecasts from [10] using 10 3 LVK events.  \n<!-- image -->  \nhave up to 171 events per year during O4. However, they use a signal-to-noise ratio (SNR) cut that is less stringent than what we use in the simulations here (which is similar to [41]), so that even if more events are expected in their case compared to [41], they are likely to have worse localization than those we consider here, so that ∼ 120 events for O4 is still a reasonable assumption for this work.  \nAs expected, constraints from runs assuming a larger value of λ tend to yield more precise measurements of the Hubble constant compared to lower λ cases, as a larger fraction of the events considered contain the 'signal' AGN flare, as opposed to background flares only, hence providing a meaningful contribution to the cosmology constraints. Also, as expected, the cases with a lower number density of AGNs often provide more stringent constraints than the case with larger number density. Perhaps more clearly, assuming a flat Ω m prior between 0 . 25 and 0 . 35 values provides more stringent constraints on H 0 . In both cases, lowering n AGN reduces the overall background count, typically resulting in better constraints on cosmological parameters, however the improvement is often within the scatter across realizations with the higher AGN density, and appears less prominent for larger values of λ . From Figure 2 it is clear that if λ ∼ 0 . 1 or greater we may reach an H 0 precision of ∼ 5% , similar to that expected from a dark siren analysis using ∼ 250 binary neutron star mergers [57], with ∼ 160 GW events.  \nIn Figure 3, we explore the more pessimistic scenarios considering values of λ of 5% , 2 . 5% and 1 . 25% . With fewer flares truly associated to BBH events compared to our fiducial case, and therefore with a larger contamination of background events, the precision increases to ∼ 6% -15% for λ = 1 . 25 -2 . 5% . This level of precision is comparable to the one available bright siren measurement [6]. It is fair to say that the BNS rates and the subsample of detectable bright sirens are highly uncertain and likely on the low end of previous estimates. Considering that for O4 at design sensitivity [41] we expect a median of 10 events with a lower limit of ∼ 0 detections and the lack of NS mergers in the first part of O4, even the scenario of λ ∼ 1% has the potential to improve upon constraints from NS mergers with only ∼ 2 ( ∼ 5 in O5) flares coming from BBHs out of 200 events at LVK design sensitivity ( ∼ 500 in O5).  \nNext, we show posterior distributions for H 0 and the dark energy equation of state. The top panel of Figure 4 presents the posteriors on H 0 with two choices of prior on Ω m after one year of A+ run, the first (dashed lines) is a flat prior Ω m ∈ [0 . 0 , 1 . 0] , the second is restricted to [0 . 25 , 0 . 35] . With this choice, we find that a ∼ 3% precision on H 0 for λ = 0 . 2 after 160 O4 events, or 2 years of O4. For comparison, we get ∼ 4 . 5% precision with the broader Ω m prior. For 500 BBH events with A+ sensitivity, we find that with the more restrictive Ω m , a 3% precision is also possible for λ = 0 . 1 .  \nThe bottom panel of Figure 4 shows some expected posteriors on the dark energy equation of state w , within a w CDM scenario. The results are for a one year run at A+ sensitivity. The w CDM model used has the same wide flat priors on Ω m ∈ [0 . 0 , 1 . 0] and a flat prior on w ∈ [ -2 . 5 , 0 . 0] . While the posteriors are broad, it is interesting that this method provides some constraining power in the Ω m -w plane. An example of a corner plot of all the parameters considered in the w CDM scenario is  \nshown in Figure 5. A banana contour typical of standard candles is seen in the Ω m -w plane, and a degeneracy is also, as expected, present between these parameters and the Hubble constant.  \nFIG. 5. Posteriors for one example run in a w CDM model with the fraction of BBH inducing flares λ = 0 . 2 , after 1 year of LVK A+ (or 500 events) considering the AGN density n AGN = 10 -4 . 75 Mpc -3 , and the fraction of AGNs flaring α flare = 10 -4 . This result also considers a flat prior over Ω m ∈ [0 . 25 , 0 . 35] .  \n<!-- image -->  \nIn Figure 6 we show the constraints we recover on λ . While these are not the main result of this manuscript, it is interesting to note how ≲ 10 -30% precision measurements are possible for all cases explored here. These constraints will have important consequences in understanding the interplay between different BBH formation channels.", 'B. AGN catalog case': 'We also apply our formalism to the case in which no follow-up is performed, therefore only an AGN catalog without flare information is used. The recovered precision on H 0 and λ AGN after 1 year ( 80 events) of LVK considering all AGNs in the GW area for n AGN = 10 -4 . 75 Mpc -3 and λ AGN in the range of [0 . 3 , 0 . 8] , is summarized in Table I. In this scenario λ AGN , differently from λ , represents the fraction of all BBHs hosted in AGNs, and thus λ AGN ≥ λ . The range of λ AGN was chosen to follow the constraints on the BBH AGN channel from [28]. We select the top 50% events in terms of comoving volume sky map localization, as that is a measure of our back-  \nTABLE I. Results in a scenario in which we apply the formalism to all AGNs in the field of a given GW event, i.e. not requiring a flare or follow-up campaigns. We present two cases, one in which the BBHs hosted in AGNs are randomly picked, and one in which the only occur in higher mass BBHs. For each case, the two rows show the H 0 precision (as median and 68% interval from the various runs made) and the recovered median fraction of BBH mergers in AGN disks, λ AGN , for 80 LVK events as a function of the input value of λ AGN (given by the different columns).  \nHigh-mass  \nσ  \n/H  \n0  \nλ  \n0  \n(%) 26  \n.  \n5  \nAGN  \n0  \n.  \n35  \n+2  \n.  \n0  \n-  \n2  \n.  \n0  \n+0  \n.  \n12  \n-  \n0  \n.  \n03  \n21  \n.  \n2  \n0  \n.  \n53  \n+5  \n.  \n9  \n-  \n5  \n.  \n3  \n+0  \n.  \n37  \n-  \n0  \n.  \n07  \n18  \n.  \n2  \n0  \n.  \n56  \n+1  \n.  \n2  \n-  \n5  \n.  \n7  \n+0  \n.  \n40  \n-  \n0  \n.  \n11  \n15  \n.  \n4  \n0  \n.  \n51  \n+5  \n.  \n5  \n-  \n4  \n.  \n0  \n+0  \n.  \n38  \n-  \n0  \n.  \n06  \n11  \n.  \n25  \n0  \n.  \n48  \n+8  \n.  \n8  \n-  \n1  \n.  \n7  \n+0  \n.  \n51  \n-  \n0  \n.  \n06  \nground contamination, and only choose the events that will provide most of the constraining power while reducing the computing time needed (which is much larger in the non-flare case compared to the flare case, given the larger contamination). Although the results show a worse precision level of H 0 compared to the flare case for a specific value of λ , as expected, they still result in an interesting precision at the ∼ 5 -10% level. Similarly to the AGN flares case, we consider the case where only BBHs on the high-mass end of the merging BH mass distribution merge in AGN disks. In that case, we find that the precision is significantly worsened, going from 5 -10% precision to ∼ 11 -27% . This is expected as the most massive mergers are more likely to be at larger distances and encompass larger localization volumes. Compared to the flares case, the contamination from background AGNs plays a more important role given that there is a difference of several orders of magnitudes, hence the worsening due to the larger volumes is more pronounced in the AGN catalog case.', 'IV. DISCUSSION': 'In this section, we discuss a few important aspects of the standard siren method discussed here. First we consider the GW selection effects, then we discuss the case of BBHs in AGNs preferentially occupying the most massive section of the BBH mass spectrum as well as caveats of the method, and finally move on to comparing to other standard siren methods. At last, we discuss some specific aspects of the case in which only AGN catalogs are avaialble instead of AGN flares.', 'A. Selection effects': 'For the flaring AGN case, we use the top 95% events based on their comoving volume localization. This cut avoids all events with localization of thousands of squared degrees. This choice is reasonable for two reasons. First,  \nthe largest localization volumes have the largest number of background contaminants, hence will provide the least constraining power to the parameters of interest. Second, given their large typical sky localization, they are unlikely to prompt follow-up observations, or be serendipitously observed by sky survey over the majority of probability of their sky localization.  \nFor the case where we take into account all AGNs, i.e. without considering the flares, we impose a more stringent cut on the localization volume. We find this to be important because in the current formalism we are assuming a certain number density of AGNs, and then assume that same value in our posterior computation. In other words, we perfectly know the expectation value of the number density of AGNs (typically down to some luminosity). As a result, perfectly knowing the number of expected AGNs within a certain comoving volume carries information about the cosmology even without having any AGN associated to the GW events (i.e. from the background AGNs only, even for input λ = 0 , effectively without a standard siren measurement). This effect has a growing impact as the background density increases, and only becomes relevant for the static AGN catalog case, as the numbers we consider for the flares are too sparse. Ideally, one would leave the AGN number density free to vary in this case, with some reasonable prior to overcome this issue and present a more realistic scenario. We test this possibility for λ AGN = 0 , and find, as expected, that it did not carry over the extra volume constraining power in the cosmological constraints. However, since it is computationally challenging to run on a large set of cases (several runs for different λ AGN values), we only constrain ourselves to the lowest background cases (i.e. lowest comoving volumes) where we have tested that the λ AGN = 0 case does not provide cosmological parameter constraints that are more constraining than the priors.  \nA systematic bias on H 0 may arise from selection effects since a specific sample of detected GW events is considered. First, to define an event detection, e.g. an SNR cut is applied. GW selection effects from this type of cuts affect a standard siren H 0 posterior approximately by a factor of ∝ H -3 0 [58], or, more generically if one considers  \nH  \n<!-- image -->  \nFIG. 6. Precision on the fraction of BBHs inducing a AGN flare, λ , as function of events in LVK (similar to the expected O4 sensitivity) and LVK A+ (similar to O5) for different true values λ . The results in the top (bottom) panel assume a flat prior over Ω m within [0 , 1 . 0] ( [0 . 25 , 0 . 35] ). The number density of AGN is n AGN = 10 -4 . 75 Mpc -3 for the filled symbols and 10 -4 . 0 Mpc -3 for the empty symbols.  \n<!-- image -->  \nother cosmological parameters, by a volume term. Because H 0 is the cosmological parameter we recover more precisely, we focus on this parameter. Since this study is focused on the precision, we disregard the influence of this selection bias in our analysis. In addition, as our formalism differs from traditional dark siren methods, our selection function includes a cosmology dependence in the µ i term in Eq. 3. A study of the impact of these selection effects on future, more precise measurements will be presented in future work. It is also worth noting that a dark standard siren H 0 posterior, and specifically the selection effects term, depend on the black hole mass distribution and merger rate redshift evolution, hence making assumptions about these (e.g. fixing the shape and parameters of the mass distribution) when computing the posterior, may introduce biases in the recovered H 0 [14].  \nOn the other hand this effect diminishes as the redshift information from the EM observations (galaxy catalog for the dark siren case, AGN observations for this work) become increasingly informative [14]. In the cases considered here, we are in a regime where the EM observations from AGNs are informative, hence it is possible that this effect does not have a significant impact given the level of statistical precision we obtain. Future work should focus on the impact of the mass distribution assumptions on the BBH-AGN formalism used here, as well as considering how to use AGN hosts in concert with spectral siren analyses as currently proposed with a generic host galaxy population [59, 60].  \nThe second selection criteria we apply is on the comoving volume. We make several tests on LVK 1 year simulations considering limited volumes (top 90% , 95% , 99% events based on volume), and without any preselection for the case of the lowest considered value of λ = 0 . 05 and highest number of background AGNs scenario with α = 10 -3 and n AGN = 10 -4 . 5 Mpc -3 and we find no significant change in the forecasts. More broadly, in our entire set of simulations we do not observe biases in the recovered posteriors when comparing to the input values (e.g. Figure 4, 5), and hence we claim that the impact of the selection effects should be smaller than the level of precision found here. On the other hand, the GW selection effects are not expected to have an impact on the recovery of λ , because, unlike the cosmology parameters, the aforementioned selection effects would not, in general, have a dependence on this parameter.  \nIn this work, we have not taken into account EM selection effects, although [39] show how these can be taken into account in the present formalism. Unlike the GW selection effects, the EM selection effects are expected to have a more significant impact on the recovery of λ rather than on the cosmology: failing to detect fainter AGNs or flares and not accounting for it, may result in recovering values of λ that are biased low compared to the truth value, but it is not expected to affect the recovered cosmology precision as the non-detection of a flare will not have constraining power on the cosmology. Since most of the events considered are at redshift z ≲ 0 . 6 (see Fig. 1), it is also reasonable to expect that current and upcoming EM facilities (e.g. Euclid, Rubin Observatory, DESI) will provide complete catalogs of AGNs at these redshifts.', 'B. BBH in AGNs as high mass mergers': 'The findings of [29] suggest that the AGN channel is more likely to produce a higher fraction of the higher mass BBHs observed in the black hole mass distribution, than those produced below the 30 -35 M ⊙ peak [56, 61, 62]. This population of mergers is typically expected to occur at a greater distance and have less precise localization than the average detected BBH. On the other hand, the results shown in Figure 2 assume that  \nFIG. 7. Hubble constant precision as a function of number of events from LVK design and LVK A+ for a scenario in which only BBH mergers with higher primary masses produce AGN flares. The upper dot-dashed horizontal line presents the expected H 0 precision from 250 Dark Siren using BNSs from [57], considering a 50% complete galaxy catalog and unweighted galaxies. The lower dot-dashed line presents the forecast from [58] after ∼ 50 BNS bright sirens. The shaded region shows the uncertainties over 100 realizations expected for a golden dark siren where we can identify the host from HLV+ with spectroscopic redshift. The number density AGN is n AGN = 10 -4 . 75 Mpc -3 for the filled symbols and 10 -4 Mpc -3 for the empty symbols. These results consider a flat prior on Ω m ∈ [0 . 25 , 0 . 35] .  \n<!-- image -->  \nthe BBHs inducing AGN flares are a random subsample of the BBH population. To take this into account, we conduct two tests. First, we assume a scenario in which BBH mergers only occur in AGN disks if the primary mass m 1 is beyond the ∼ 30 -35 M ⊙ peak observed in the mass distribution. We assume m 1 > 35 M ⊙ for the BBHs in AGNs for simplicity, as it is reasonable to believe that the peak may originate through a formation channel different from the AGN one, and this cut is justified considering that [62] finds a slightly lower value of the peak mass than previously found for GWTC-3. We present the results for this case in Figure 7. We find that the precision worsens, as expected, compared to our previous findings, although only slightly when comparing results for the same values of λ . This is due to the fact that the number of events occurring in AGNs, hence λ , is diminished (even if we assume that 100% BBHs with m 1 > 35 M ⊙ produce a flare, we cannot reach the λ = 0.2 or 0.3 cases in Figure 2), and because these massive mergers will typically have worse localizations and distance precision. Comparing e.g. the λ = 0 . 1 case of Figure 2 with that of Figure 7, we find slightly larger scatter in the H 0 precision in the latter case, as even a single nearby event detection could be very loud and provide better distance measurements than lower mass objects at a similar distance. For the second test, we  \nalso consider lower mass BBHs. If we consider that the higher mass BH in the mass distributions are formed hierarchically in the AGN disk, then it is also reasonable to assume that for each BBH with m 1 > 35 M ⊙ , there is at least another BBH of first generation black holes below that value that formed the second generation BH in the m 1 > 35 M ⊙ binary, so we also sample for each high mass BBH in AGN a low mass BBH in AGN with primary below the peak. The results in this case become extremely similar to these in Figure 2 for the same values of λ , so we do not show them in the figure. We conclude that a 2-10% precision in H 0 is also possible with higher mass BBHs merging in AGNs with a A+ sensitivity, assuming that λ > 1% .', 'C. Caveats': 'A caveat of the proposed method is that its performance is closely related to the fraction of BBH merging in AGN disks. GW population studies [63] show that multiple BBH formation channels are likely at play, with dynamical formation being a viable scenario, but the contribution of each channel is yet to be well constrained. It is also worth noticing that our method, at least for the AGN flare case, is highly dependent upon the detectability of EM counterparts of BBHs inducing AGN flares, which is a subfraction of the BBHs hosted in AGN disks. The detectability is closely related to the interaction of BBHs and the AGN disks proprieties and is being currently investigated. This hypothesis has theoretical modelling support [31-35], and promising candidates from these channels have been reported [36, 37]. For instance, [37] identified several new candidate EM counterparts to BBH events during the LVK O3 campaign. According to the authors, 7 EM counterparts were statistically associated to BBH events (with a coincidence probability of the flares with all the associated BBH of p = 1 . 90 × 10 -3 ) from the total set of 83 events considered on the third observing run (O3). Their findings suggest that the fraction of detectable flares from BBHs is around ∼ 0 . 07 in O3 within the ZTF explored volume up to magnitude g < 20 . 5 , considering ZTF observes roughly half of the sky.  \nAnother caveat for the forecast presented here is the assumption that some fraction of the observed BBH population occurs in AGN, where there is no correlation between this specific formation channel and the BBH properties. In reality, it is possible that the AGN formation channel dominates the higher mass end of the BH mass distribution, thanks to its ability to retain hierarchical mergers.  \nIt is important to mention that the forecasts presented here (as well as in most of the other standard siren forecast works quoted) assume the LVK detectors at a sensitivity (roughly the design sensitivity) that was expected until mid 2022 for the upcoming observing run O4. In reality, the Virgo detector has not joined the O4 run, and  \nit is unclear whether its expected sensitivity for O4 will be achieved during this run. KAGRA has joined the run for a limited time at a significantly lower BNS inspiral range than what assumed here (80 Mpc). These changes to the observing capabilities only became publicly clear after this work was close to finalization. Therefore, the forecasts presented here are optimistic if considered for the upcoming O4 run, but they are still useful for future runs.', 'D. Comparison with other standard siren methods': 'Next, we consider how competitive and complementary our method is when compared to several other standard siren constraints expected from current generation GW detectors.  \nFirst we consider the constraints shown for reference by the horizontal lines in Figure 2. We stress that this is not a realistic comparison in terms of upcoming results since the analysis of [57] is based on BNS simulations that assume an O2-like sensitivity (hence, the same event detected in e.g. O4 should have an improved signal-tonoise ratio, hence localization, compared to those simulations on average). Moreover, 250 BNSs are unlikely to be detected during O4 (and potentially not even in 1 year of O5). Note that [58] provides less optimistic forecasts from dark standard sirens with BBH mergers, potentially only reaching a 10% precision with all detectors considered here at design sensitivity, while a BNS dark siren analysis is expected to reach that 5% within such scenario. On the other hand, [19] show how some extremely well-localized dark sirens could provide a ∼ 2 -8% precision (shown as shaded region in Figure 2). It is also worth mentioning that a precision of ∼ 2 . 5 -3% is feasible within about a year of A+ if λ is on the high-end ( λ ∼ 0 . 2 -0 . 3 ) of the range we consider. [57] consider BNS mergers dark sirens which are typically closer than BBHs and less sensitive to Ω m . For a spectroscopic galaxy catalog which is 50% complete, they find a 5% precision on H 0 , which is a constraint that our method can also achieve within the next ∼ 5 years as long as λ > 0 . 1 . The study in [58] relies primarily on bright standard sirens, and finds that to reach the 2% precision level would take around 50 bright sirens. Considering the median detections expected in [51] scenarios, these figures are likely to be reached during O5. On the other hand, considering that only one BNS EM counterpart has been detected so far and that follow-up campaigns might not succeed on detecting all EM counterparts, due to limited telescope time or more edge-on or high mass systems dominated by a fainter and/or red component [52], the number of events with confirmed counterparts may be significantly less than the number of BNS or NSBH detections. The addition of another sample of candidate EM counterparts and potential standard sirens is therefore interesting to ensure that GW cosmology is able to weigh in the Hubble constant tension before the 2030s.  \nCompared to other bright standard siren classes, namely BNS and NSBH, the observation of flares might not require disruptive Target-of-Opportunity (ToO) imaging or even ToO spectra for the present study. First, there is no need to require a direct association of a given flare to a specific GW event, although ideally spectroscopic confirmation of e.g. off-center events through spectroscopy of asymmetric broad lines [32] would be a smoking gun for association (hence becoming a bright standard siren, if the uniqueness of the counterpart can be established). Additionally, the flare time scales are expected to be of the orders of weeks to months for higher mass AGNs, as opposed to less than 1 week for kilonovae. Hence, AGN flare counterpart searches would rather require all sky monitoring such as that carried out by ZTF or Rubin Legacy Survey of Space and Time (LSST), and a small spectroscopic follow-up campaign to acquire redshifts of flares where not available. For instance, in our high background contaminant scenario, our simulations show that ∼ 95% of events with n AGN = 10 -4 . 5 Mpc -3 and α flare = 10 -3 would have less than ∼ 7000 AGNs in the GW area and therefore around 7 background flares potentially associated to an event.  \n[19] requires golden events where one can identify the exact BBH host galaxy, assuming 100% duty cycle over 2 years of LVK at A+ sensitivity, to reach a ∼ 2% precision level on H 0 . Under these circumstances, dark sirens act as bright sirens, and are therefore extremely powerful. In these special cases, no active follow-up is required for the purpose of dark standard siren cosmology, although complete imaging and spectroscopy of the high probability region will be fundamental. On the other hand, discovering an AGN flare in coincidence with such golden events would constitute a smoking gun for the GW-AGN association and a higher confidence standard siren measurement, given the low probability of chance coincidence of flares over the small volumes considered.  \n[10] present a method that constrains the Hubble parameter using GW BBHs alone by taking advantage of a specific feature in the mass distribution of black holes. Their method is most constraining at H ( z = 0 . 8) with a set of ∼ 10 3 to 5 × 10 3 events from the LVK network and can be translated to a precision in H 0 to ∼ 6% ( ∼ 10 3 events) and ∼ 3% ( ∼ 5 × 10 3 ) level. The method presented in this work is competitive, considering a free Ω m ∈ [0 . 0 , 1 . 0] we use 240 events for LVK to reach ∼ 5% level for λ = 0 . 1 . However, the method of [10] has the advantage of not requiring any GW follow-up. The expectation of the H ( z ) analysis is that it will enable constraints on w to a 19% ( ∼ 10 3 events) and 12% ( ∼ 5 × 10 3 events) precision after 1 and 5 years of LVK A+ (O5) running at design sensitivity, respectively, after imposing a 1% prior on H 0 and the Planck 2018 constraint on Ω m . Additionally, [11] find slightly less constraining forecasts on H ( z ) with an updated BBH population. Although we present wider uncertainties on w our method is still promising compared with these proposed analyses, and perhaps more importantly, it is complementary. While  \nmost of the constraining power for the spectral sirens typically comes from the redshift and mass regions where the highest statistics of detections are available, specifically for 2G detectors, the lower bound of the pair-instability mass gap [11], or potentially the ∼ 30 M ⊙ peak, and distances close to the detection horizon, in our case the constraints are more powerful at lower distances where better localizations are more likely. It will eventually be interesting to combine the two methods.', 'E. AGN catalog case': 'At last, we discuss the case of follow-up for the AGN catalog case, where the flares are not taken into account. Considering only the best localized events will limit the number of AGNs in the field for which spectroscopic redshifts are required to ∼ 50 -70 at maximum per event, and thus feasible within the reach of many current spectroscopic facilities for a dedicated campaign. Such events are also the most relevant events since they have a lower number of background contaminants and lower uncertainties in d L . The usage of the present methodology with the full set of AGNs in the field can be interpreted as a type of informed dark standard siren method, which has the advantage of not requiring full galaxy catalogs, but the disadvantage of not providing constraining power for events that do not occur in AGNs. Our forecasts show that a 5-10% precision on H 0 could be achieved within 1 year of LVK at design sensitivity in the optimistic case that 30-80% of BBHs occur in AGNs. In the future, it will be interesting to expand on this promising avenue by taking into account AGN clustering in the simulations, and by further speeding up the sampling to produce full ensembles of simulated scenarios.  \nWe also briefly discuss the role of completeness of AGN catalogs. For what concerns the AGN flares follow-up, although complete AGN catalogs would facilitate galaxytargeted searches, they are not necessary for the method. For wide field searches, such as those that can be carried out with ZTF and LSST, it is reasonable to search for nuclear transients regardless of existing catalogs. Since nuclear flares of the type and shape that we would search for are typically rare, if some of the flares do not have an AGN classification, one could take additional observations to clarify their nature. Similarly, one may consider to observe the localization of the GW event if a BBH shows features that may hint to an AGN formation, to produce an AGN catalog down to the desired luminosities. With the upcoming completion of large spectroscopic surveys such as DESI [64], and the start of new surveys [65], it is possible that dedicated follow-up will not be needed as complete catalogs of Type I AGNs will be available. For both the flare and the AGN catalog case, the formalism allows us to specify the selection func-  \ntion of AGNs and flares (see Eq. 3). For both cases, it will be interesting in the future to quantify which AGNs luminosities are likely to contribute significantly to the BBH merger rates, and apply the selection function for those (efforts in that direction already exist, e.g. [66]).', 'V. CONCLUSIONS': 'In this study, we present forecasts using a novel and competitive method to obtain cosmological constraints from multi-messenger observations of gravitational wave events. Our work is build upon the formalism proposed in [39], where the perspective of using AGN flares to constrain cosmology is pointed out. From a set of simulations representing a possible sample of detections from the upcoming LVK observing runs, we obtain constraints for cosmological parameters within the Λ CDM and w CDM models using multiple realizations.  \nWe find interesting constraints on the Hubble constant are possible, down to 3 -6% precision for O4 and O5-like runs, assuming that at least ∼ 10% of BBHs occurs in AGN disks. We reach more competitive constraints on H 0 , as opposed to Ω m or w , because most of the constraining power comes from the events localized within a smaller volume, typically the most nearby, at z < 0 . 2 . The constraints we find will of course only be competitive if a significant fraction of LVK BBH mergers occur in AGNs, which is yet to be confirmed. Assuming that the AGN origin of BBHs can be confirmed, future work on this subject should explore in more detail how selection effects affect standard siren constraints from the proposed method.', 'ACKNOWLEDGMENTS': "Antonella Palmese thanks Maya Fishbach, Saavik Ford, Matthew Graham, Barry McKernan, Saul Perlmutter for useful discussion, Leo Singer for help with the injections and sky map software and for making BAYESTAR public, and the LIGO collaboration for providing public access tools. Antonella Palmese acknowledges support for this work was provided by NASA through the NASA Hubble Fellowship grant HST-HF251488.001-A awarded by the Space Telescope Science Institute, which is operated by Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. Clecio Bom 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). 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.  \n- [1] B. F. Schutz, Determining the Hubble constant from\n- gravitational wave observations, Nature 323 , 310 (1986). [2] Planck Collaboration, N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, S. Basak, R. Battye, K. Benabed, J. P. Bernard, M. Bersanelli, P. Bielewicz, J. J. Bock, J. R. Bond, J. Borrill, F. R. Bouchet, F. Boulanger, M. Bucher, C. Burigana, R. C. Butler, E. Calabrese, J. F. Cardoso, J. Carron, A. Challinor, H. C. Chiang, J. Chluba, L. P. L. Colombo, C. Combet, D. Contreras, B. P. Crill, F. Cuttaia, P. de Bernardis, G. de Zotti, J. Delabrouille, J. M. Delouis, E. Di Valentino, J. M. Diego, O. Doré, M. Douspis, A. Ducout, X. Dupac, S. Dusini, G. Efstathiou, F. Elsner, T. A. Enßlin, H. K. Eriksen, Y. Fantaye, M. Farhang, J. Fergusson, R. Fernandez-Cobos, F. Finelli, F. Forastieri, M. Frailis, A. A. Fraisse, E. Franceschi, A. Frolov, S. Galeotta, S. Galli, K. Ganga, R. T. Génova-Santos, M. Gerbino, T. Ghosh, J. González-Nuevo, K. M. Górski, S. Gratton, A. Gruppuso, J. E. Gudmundsson, J. Hamann, W. Handley, F. K. Hansen, D. Herranz, S. R. Hildebrandt, E. Hivon, Z. Huang, A. H. Jaffe, W. C. Jones, A. Karakci, E. Keihänen, R. Keskitalo, K. Kiiveri, J. Kim, T. S. Kisner, L. Knox, N. Krachmalnicoff, M. Kunz, H. Kurki-Suonio, G. Lagache, J. M. Lamarre, A. Lasenby, M. Lattanzi, C. R. Lawrence, M. Le Jeune, P. Lemos, J. Lesgourgues, F. Levrier, A. Lewis, M. Liguori, P. B. Lilje, M. Lilley, V. Lindholm, M. López-Caniego, P. M. Lubin, Y. Z. Ma, J. F. Macías-Pérez, G. Maggio, D. Maino, N. Mandolesi, A. Mangilli, A. Marcos-Caballero, M. Maris, P. G. Martin, M. Martinelli, E. Martínez-González, S. Matarrese, N. Mauri, J. D. McEwen, P. R. Meinhold, A. Melchiorri, A. Mennella, M. Migliaccio, M. Millea, S. Mitra, M. A. Miville-Deschênes, D. Molinari, L. Montier, G. Morgante, A. Moss, P. Natoli, H. U. Nørgaard-Nielsen, L. Pagano, D. Paoletti, B. Partridge, G. Patanchon, H. V. Peiris, F. Perrotta, V. Pettorino, F. Piacentini, L. Polastri, G. Polenta, J. L. Puget, J. P. Rachen, M. Reinecke, M. Remazeilles, A. Renzi, G. Rocha, C. Rosset, G. Roudier, J. A. Rubiño-Martín, B. Ruiz-Granados, L. Salvati, M. Sandri, M. Savelainen, D. Scott, E. P. S. Shellard, C. Sirignano, G. Sirri, L. D. Spencer, R. Sunyaev, A. S. Suur-Uski, J. A. Tauber, D. Tavagnacco, M. Tenti, L. Toffolatti, M. Tomasi, T. Trombetti, L. Valenziano, J. Valiviita, B. Van Tent, L. Vibert, P. Vielva, F. Villa, N. Vittorio, B. D. Wandelt, I. K. Wehus, M. White, S. D. M. White, A. Zacchei, and A. Zonca, Planck 2018 results. VI. Cosmological parameters, A&A 641 , A6 (2020), arXiv:1807.06209 [astro-ph.CO].\n- [3] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Cosmic distances calibrated to 1% precision with gaia edr3 parallaxes and hubble space telescope photometry of 75 milky way cepheids confirm tension with λ cdm, The Astrophysical Journal 908 , L6 (2021).\n- [4] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119 , 161101 (2017).  \n- [5] LIGO Scientific Collaboration, Virgo Collaboration, F. GBM, INTEGRAL, IceCube Collaboration, AstroSat Cadmium Zinc Telluride Imager Team, IPN Collaboration, The Insight-Hxmt Collaboration, ANTARES Collaboration, The Swift Collaboration, AGILE Team, The 1M2H Team, The Dark Energy Camera GW-EM Collaboration, the DES Collaboration, The DLT40 Collaboration, GRAWITA, :, GRAvitational Wave Inaf TeAm, The Fermi Large Area Telescope Collaboration, ATCA, :, A. Telescope Compact Array, ASKAP, :, A. SKA Pathfinder, Las Cumbres Observatory Group, OzGrav, DWF, AST3, CAASTRO Collaborations, The VINROUGE Collaboration, MASTER Collaboration, JGEM, GROWTH, JAGWAR, C. NRAO, TTU-NRAO, NuSTAR Collaborations, Pan-STARRS, The MAXI Team, T. Consortium, KU Collaboration, N. Optical Telescope, ePESSTO, GROND, T. Tech University, SALT Group, TOROS, :, Transient Robotic Observatory of the South Collaboration, The BOOTES Collaboration, MWA, :, M. Widefield Array, The CALET Collaboration, IKI-GW Follow-up Collaboration, H. E. S. S. Collaboration, LOFAR Collaboration, LWA, :, L. Wavelength Array, HAWC Collaboration, The Pierre Auger Collaboration, ALMA Collaboration, Euro VLBI Team, Pi of the Sky Collaboration, The Chandra Team at McGill University, DFN, :, D. Fireball Network, ATLAS, H. Time Resolution Universe Survey, RIMAS, RATIR, and S. South Africa/MeerKAT, Multi-messenger Observations of a Binary Neutron Star Merger, ApJ 848 , L12 (2017), arXiv:1710.05833 [astro-ph.HE].\n- [6] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, and et al., A gravitational-wave standard siren measurement of the Hubble constant, Nature 551 , 85 (2017), arXiv:1710.05835.\n- [7] A. Palmese, R. Kaur, A. Hajela, R. Margutti, A. McDowell, and A. MacFadyen, A standard siren measurement of the Hubble constant using GW170817 and the latest observations of the electromagnetic counterpart afterglow, arXiv e-prints , arXiv:2305.19914 (2023), arXiv:2305.19914 [astro-ph.CO].\n- [8] K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K. P. Mooley, and A. T. Deller, A Hubble constant measurement from superluminal motion of the jet in GW170817, Nature Astronomy 3 , 940 (2019), arXiv:1806.10596 [astro-ph.CO].\n- [9] S. M. Feeney, H. V. Peiris, A. R. Williamson, S. M. Nissanke, D. J. Mortlock, J. Alsing, and D. Scolnic, Prospects for Resolving the Hubble Constant Tension with Standard Sirens, Phys. Rev. Lett. 122 , 061105 (2019), arXiv:1802.03404 [astro-ph.CO].\n- [10] W. M. Farr, M. Fishbach, J. Ye, and D. E. Holz, A future percent-level measurement of the hubble expansion at redshift 0.8 with advanced ligo, ApJ 883 , L42 (2019).\n- [11] J. M. Ezquiaga and D. E. Holz, Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries, Phys. Rev. Lett. 129 , 061102 (2022), arXiv:2202.08240 [astro-ph.CO].\n- [12] M. Soares-Santos, A. Palmese, W. Hartley, J. Annis, J. Garcia-Bellido, O. Lahav, Z. Doctor, M. Fishbach, D. E. Holz, H. Lin, and et al., First Measure-  \nment of the Hubble Constant from a Dark Standard Siren using the Dark Energy Survey Galaxies and the LIGO/Virgo Binary-Black-hole Merger GW170814, ApJ 876 , L7 (2019), arXiv:1901.01540.  \n- [13] J. R. Gair, A. Ghosh, R. Gray, D. E. Holz, S. Mastrogiovanni, S. Mukherjee, A. Palmese, N. Tamanini, T. Baker, F. Beirnaert, M. Bilicki, H.-Y. Chen, G. Dálya, J. M. Ezquiaga, W. M. Farr, M. Fishbach, J. GarciaBellido, T. Ghosh, H.-Y. Huang, C. Karathanasis, K. Leyde, I. Magaña Hernandez, J. Noller, G. Pierra, P. Raffai, A. Enea Romano, M. Seglar-Arroyo, D. A. Steer, C. Turski, M. P. Vaccaro, and S. Andrés VallejoPeña, The Hitchhiker's guide to the galaxy catalog approach for gravitational wave cosmology, arXiv e-prints , arXiv:2212.08694 (2022), arXiv:2212.08694 [gr-qc].\n- [14] LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration, Constraints on the cosmic expansion history from gwtc-3 (2021), arXiv:2111.03604 [astro-ph.CO].\n- [15] A. Palmese, C. R. Bom, S. Mucesh, and W. G. Hartley, A standard siren measurement of the Hubble constant using gravitational wave events from the first three LIGO/Virgo observing runs and the DESI Legacy Survey, arXiv e-prints , arXiv:2111.06445 (2021), arXiv:2111.06445 [astro-ph.CO].\n- [16] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and M. Mancarella, Cosmology with ligo/virgo dark sirens: Hubble parameter and modified gravitational wave propagation (2021), arXiv:2101.12660 [astro-ph.CO].\n- [17] A. Palmese, J. deVicente, M. E. S. Pereira, J. Annis, W. Hartley, K. Herner, M. Soares-Santos, M. Crocce, D. Huterer, I. Magaña Hernandez, A. Garcia, J. GarciaBellido, J. Gschwend, D. E. Holz, R. Kessler, O. Lahav, R. Morgan, C. Nicolaou, C. Conselice, R. J. Foley, M. S. S. Gill, T. M. C. Abbott, M. Aguena, S. Allam, S. Avila, K. Bechtol, E. Bertin, S. Bhargava, D. Brooks, E. Buckley-Geer, D. L. Burke, M. Carrasco Kind, J. Carretero, F. J. Castand er, C. Chang, M. Costanzi, L. N. da Costa, T. M. Davis, S. Desai, H. T. Diehl, P. Doel, A. Drlica-Wagner, J. Estrada, S. Everett, A. E. Evrard, E. Fernandez, D. A. Finley, B. Flaugher, P. Fosalba, J. Frieman, E. Gaztanaga, D. W. Gerdes, D. Gruen, R. A. Gruendl, G. Gutierrez, S. R. Hinton, D. L. Hollowood, K. Honscheid, D. J. James, S. Kent, E. Krause, K. Kuehn, H. Lin, M. A. G. Maia, M. March, J. L. Marshall, P. Melchior, F. Menanteau, R. Miquel, R. L. C. Ogando, F. Paz-Chinchón, A. A. Plazas, A. Roodman, M. Sako, E. Sanchez, V. Scarpine, M. Schubnell, S. Serrano, I. Sevilla-Noarbe, J. A. Smith, M. Smith, E. Suchyta, G. Tarle, M. A. Troxel, D. L. Tucker, A. R. Walker, W. Wester, R. D. Wilkinson, J. Zuntz, and DES Collaboration, A Statistical Standard Siren Measurement of the Hubble Constant from the LIGO/Virgo Gravitational Wave Compact Object Merger GW190814 and Dark Energy Survey Galaxies, ApJ 900 , L33 (2020), arXiv:2006.14961 [astro-ph.CO].\n- [18] W. Ballard, A. Palmese, I. M. Hernandez, S. BenZvi, J. Moon, A. J. Ross, G. Rossi, J. Aguilar, S. Ahlen, R. Blum, D. Brooks, T. Claybaugh, A. de la Macorra, A. Dey, P. Doel, J. E. Forero-Romero, S. Gontcho A Gontcho, K. Honscheid, A. Kremin, M. Manera, A. Meisner, R. Miquel, J. Moustakas, F. Prada, E. Sanchez, G. Tarlé, Z. Zhou, and DESI Collaboration, A Dark Siren Measurement of the Hubble Constant with the  \n- LIGO/Virgo Gravitational Wave Event GW190412 and DESI Galaxies, Research Notes of the American Astronomical Society 7 , 250 (2023), arXiv:2311.13062 [astroph.CO].\n- [19] S. Borhanian, A. Dhani, A. Gupta, K. G. Arun, and B. S. Sathyaprakash, Dark Sirens to Resolve the HubbleLemaître Tension, arXiv e-prints , arXiv:2007.02883 (2020), arXiv:2007.02883 [astro-ph.CO].\n- [20] C. Cigarrán Díaz and S. Mukherjee, Mapping the cosmic expansion history from LIGO-Virgo-KAGRA in synergy with DESI and SPHEREx, MNRAS 511 , 2782 (2022), arXiv:2107.12787 [astro-ph.CO].\n- [21] K. Belczynski, V. Kalogera, and T. Bulik, A comprehensive study of binary compact objects as gravitational wave sources: Evolutionary channels, rates, and physical properties, The Astrophysical Journal 572 , 407 (2002).\n- [22] C. L. Rodriguez, C.-J. Haster, S. Chatterjee, V. Kalogera, and F. A. Rasio, DYNAMICAL FORMATION OF THE GW150914 BINARY BLACK HOLE, The Astrophysical Journal 824 , L8 (2016).\n- [23] B. McKernan, K. E. S. Ford, W. Lyra, and H. B. Perets, Intermediate mass black holes in AGN discs I. Production and growth, MNRAS 425 , 460 (2012), arXiv:1206.2309 [astro-ph.GA].\n- [24] Y. Yang, I. Bartos, V. Gayathri, K. E. S. Ford, Z. Haiman, S. Klimenko, B. Kocsis, S. Márka, Z. Márka, B. McKernan, and R. O'Shaughnessy, Hierarchical Black Hole Mergers in Active Galactic Nuclei, Phys. Rev. Lett. 123 , 181101 (2019), arXiv:1906.09281 [astro-ph.HE].\n- [25] C. J. Conselice, R. Bhatawdekar, A. Palmese, and W. G. Hartley, Gravitational Waves from Black Holes in Merging Ultra-Dwarf Galaxies, ApJ 890 , arXiv:1907.05361 (2020), arXiv:1907.05361 [astro-ph.GA].\n- [26] A. Palmese and C. J. Conselice, GW190521 from the Merger of Ultradwarf Galaxies, Phys. Rev. Lett. 126 , 181103 (2021), arXiv:2009.10688 [astro-ph.GA].\n- [27] Y. Yang, I. Bartos, Z. Haiman, B. Kocsis, Z. Márka, N. C. Stone, and S. Márka, AGN Disks Harden the Mass Distribution of Stellar-mass Binary Black Hole Mergers, ApJ 876 , 122 (2019), arXiv:1903.01405 [astro-ph.HE].\n- [28] K. E. S. Ford and B. McKernan, Binary black hole merger rates in AGN discs versus nuclear star clusters: loud beats quiet, MNRAS 517 , 5827 (2022), arXiv:2109.03212 [astro-ph.HE].\n- [29] V. Gayathri, Y. Yang, H. Tagawa, Z. Haiman, and I. Bartos, Black hole mergers of AGN origin in LIGO-virgo's o1-o3a observing periods, The Astrophysical Journal Letters 920 , L42 (2021).\n- [30] V. Gayathri, D. Wysocki, Y. Yang, V. Delfavero, R. O'Shaughnessy, Z. Haiman, H. Tagawa, and I. Bartos, Gravitational Wave Source Populations: Disentangling an AGN Component, ApJ 945 , L29 (2023), arXiv:2301.04187 [gr-qc].\n- [31] I. Bartos, B. Kocsis, Z. Haiman, and S. Márka, Rapid and Bright Stellar-mass Binary Black Hole Mergers in Active Galactic Nuclei, ApJ 835 , 165 (2017), arXiv:1602.03831 [astro-ph.HE].\n- [32] B. McKernan, K. E. S. Ford, I. Bartos, M. J. Graham, W. Lyra, S. Marka, Z. Marka, N. P. Ross, D. Stern, and Y. Yang, Ram-pressure Stripping of a Kicked Hill Sphere: Prompt Electromagnetic Emission from the Merger of Stellar Mass Black Holes in an AGN Accretion Disk, ApJ 884 , L50 (2019), arXiv:1907.03746 [astro-ph.HE].\n- [33] S. S. Kimura, K. Murase, and I. Bartos, Outflow Bubbles from Compact Binary Mergers Embedded in Active Galactic Nuclei: Cavity Formation and the Impact on Electromagnetic Counterparts, ApJ 916 , 111 (2021), arXiv:2103.02461 [astro-ph.HE].\n- [34] H. Tagawa, S. S. Kimura, Z. Haiman, R. Perna, and I. Bartos, Observable signatures of stellar-mass black holes in active galactic nuclei, arXiv e-prints , arXiv:2303.02172 (2023), arXiv:2303.02172 [astroph.HE].\n- [35] J. Rodríguez-Ramírez, C. Bom, B. Fraga, and R. Nemmen, Optical emission model for binary black hole merger remnants travelling through discs of active galactic nucleus, arXiv preprint arXiv:2304.10567 (2023).\n- [36] M. J. Graham, K. E. S. Ford, B. McKernan, N. P. Ross, D. Stern, K. Burdge, M. Coughlin, S. G. Djorgovski, A. J. Drake, D. Duev, M. Kasliwal, A. A. Mahabal, S. van Velzen, J. Belecki, E. C. Bellm, R. Burruss, S. B. Cenko, V. Cunningham, G. Helou, S. R. Kulkarni, F. J. Masci, T. Prince, D. Reiley, H. Rodriguez, B. Rusholme, R. M. Smith, and M. T. Soumagnac, Candidate Electromagnetic Counterpart to the Binary Black Hole Merger Gravitational-Wave Event S190521g ∗ , Phys. Rev. Lett. 124 , 251102 (2020), arXiv:2006.14122 [astro-ph.HE].\n- [37] M. J. Graham, B. McKernan, K. E. S. Ford, D. Stern, S. G. Djorgovski, M. Coughlin, K. B. Burdge, E. C. Bellm, G. Helou, A. A. Mahabal, F. J. Masci, J. Purdum, P. Rosnet, and B. Rusholme, A Light in the Dark: Searching for Electromagnetic Counterparts to Black Hole-Black Hole Mergers in LIGO/Virgo O3 with the Zwicky Transient Facility, ApJ 942 , 99 (2023), arXiv:2209.13004 [astro-ph.HE].\n- [38] I. Bartos, Z. Haiman, Z. Marka, B. D. Metzger, N. C. Stone, and S. Marka, Gravitational-wave localization alone can probe origin of stellar-mass black hole mergers, Nature Communications 8 , 10.1038/s41467-017-00851-7 (2017).\n- [39] A. Palmese, M. Fishbach, C. J. Burke, J. Annis, and X. Liu, Do ligo/virgo black hole mergers produce agn flares? the case of gw190521 and prospects for reaching a confident association, The Astrophysical Journal Letters 914 , L34 (2021).\n- [40] N. Veronesi, E. M. Rossi, S. van Velzen, and R. Buscicchio, Detectability of a spatial correlation between stellar mass black hole mergers and active galactic nuclei in the local Universe, MNRAS 514 , 2092 (2022), arXiv:2203.05907 [astro-ph.HE].\n- [41] B. P. Abbott, R. Abbott, T. D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, V. B. Adya, C. Affeldt, M. Agathos, K. Agatsuma, N. Aggarwal, O. D. Aguiar, L. Aiello, A. Ain, P. Ajith, T. Akutsu, G. Allen, A. Allocca, M. A. Aloy, P. A. Altin, A. Amato, A. Ananyeva, S. B. Anderson, W. G. Anderson, M. Ando, S. V. Angelova, S. Antier, S. Appert, K. Arai, K. Arai, Y. Arai, S. Araki, A. Araya, M. C. Araya, J. S. Areeda, M. Arene, N. Aritomi, N. Arnaud, K. G. Arun, S. Ascenzi, G. Ashton, Y. Aso, S. M. Aston, P. Astone, F. Aubin, P. Aufmuth, K. AultONeal, C. Austin, V. Avendano, A. Avila-Alvarez, S. Babak, P. Bacon, F. Badaracco, M. K. M. Bader, S. W. Bae, Y. B. Bae, L. Baiotti, and Others, Prospects for observing and localizing gravitational-wave transients with advanced ligo, advanced virgo and kagra, Living Reviews in Relativity 23 , 10.1007/s41114-020-00026-9 (2020).  \n- [42] Y. Yang, V. Gayathri, L. M. B. Alves, A. G. Sullivan, Z. Marka, S. Marka, and I. Bartos, Determining the Hubble Constant with AGN-assisted Black Hole Mergers, arXiv e-prints , arXiv:2009.13739 (2020), arXiv:2009.13739 [astro-ph.HE].\n- [43] S. Mukherjee, A. Ghosh, M. J. Graham, C. Karathanasis, M. M. Kasliwal, I. Magaña Hernandez, S. M. Nissanke, A. Silvestri, and B. D. Wandelt, First measurement of the Hubble parameter from bright binary black hole GW190521, arXiv e-prints , arXiv:2009.14199 (2020), arXiv:2009.14199 [astro-ph.CO].\n- [44] V. Gayathri, J. Healy, J. Lange, B. O'Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. O. Lousto, and R. O'Shaughnessy, Measuring the Hubble Constant with GW190521 as an Eccentric black hole Merger and Its Potential Electromagnetic Counterpart, ApJ 908 , L34 (2021).\n- [45] H.-Y. Chen, C.-J. Haster, S. Vitale, W. M. Farr, and M. Isi, A standard siren cosmological measurement from the potential gw190521 electromagnetic counterpart ztf19abanrhr, Monthly Notices of the Royal Astronomical Society 513 , 2152-2157 (2022).\n- [46] J. E. Greene and L. C. Ho, The mass function of active black holes in the local universe, The Astrophysical Journal 667 , 131 (2007).\n- [47] J. E. Greene and L. C. Ho, Erratum: 'the mass function of active black holes in the local universe' (2007, apj, 667, 131), The Astrophysical Journal 704 , 1743 (2009).\n- [48] I. Bartos, Z. Haiman, Z. Marka, B. D. Metzger, N. C. Stone, and S. Marka, Gravitational-wave localization alone can probe origin of stellar-mass black hole mergers, Nature Communications 8 , 831 (2017), arXiv:1701.02328 [astro-ph.HE].\n- [49] E. C. Kool, T. M. Reynolds, S. Mattila, E. Kankare, M. A. Pérez-Torres, A. Efstathiou, S. Ryder, C. RomeroCañizales, W. Lu, T. Heikkilä, G. E. Anderson, M. Berton, J. Bright, G. Cannizzaro, D. Eappachen, M. Fraser, M. Gromadzki, P. G. Jonker, H. Kuncarayakti, P. Lundqvist, K. Maeda, R. M. McDermid, A. M. Medling, S. Moran, A. Reguitti, M. Shahbandeh, S. Tsygankov, V. U, and T. Wevers, AT 2017gbl: a dust obscured TDE candidate in a luminous infrared galaxy, Monthly Notices of the Royal Astronomical Society 498 , 2167 (2020), https://academic.oup.com/mnras/articlepdf/498/2/2167/33776981/staa2351.pdf.\n- [50] P. F. Hopkins, G. T. Richards, and L. Hernquist, An observational determination of the bolometric quasar luminosity function, The Astrophysical Journal 654 , 731 (2007).\n- [51] P. Petrov, L. P. Singer, M. W. Coughlin, V. Kumar, M. Almualla, S. Anand, M. Bulla, T. Dietrich, F. Foucart, and N. Guessoum, Data-driven expectations for electromagnetic counterpart searches based on ligo/virgo public alerts, The Astrophysical Journal 924 , 54 (2022).\n- [52] C. R. Bom, J. Annis, A. Garcia, A. Palmese, N. Sherman, M. Soares-Santos, L. Santana-Silva, R. Morgan, K. Bechtol, T. Davis, et al. , Designing an optimal kilonova search using decam for gravitational wave events, arXiv preprint arXiv:2302.04878 (2023).\n- [53] L. P. Singer, H.-Y. Chen, D. E. Holz, W. M. Farr, L. R. Price, V. Raymond, S. B. Cenko, N. Gehrels, J. Cannizzo, M. M. Kasliwal, S. Nissanke, M. Coughlin, B. Farr, A. L. Urban, S. Vitale, J. Veitch, P. Graff, C. P. L. Berry, S. Mohapatra, and I. Mandel, GOING\n- THE DISTANCE: MAPPING HOST GALAXIES OF LIGO AND VIRGO SOURCES IN THREE DIMENSIONS USING LOCAL COSMOGRAPHY AND TARGETED FOLLOW-UP, The Astrophysical Journal 829 , L15 (2016).\n- [54] L. P. Singer and L. R. Price, Rapid bayesian position reconstruction for gravitational-wave transients, Phys. Rev. D 93 , 024013 (2016).\n- [55] LIGO Scientific Collaboration, LIGO Algorithm Library - LALSuite, free software (GPL) (2018).\n- [56] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13 , 011048 (2023), arXiv:2111.03634 [astro-ph.HE].\n- [57] R. Gray, I. M. Hernandez, H. Qi, A. Sur, P. R. Brady, H.Y. Chen, W. M. Farr, M. Fishbach, J. R. Gair, A. Ghosh, D. E. Holz, S. Mastrogiovanni, C. Messenger, D. A. Steer, and J. Veitch, Cosmological inference using gravitational wave standard sirens: A mock data analysis, Phys. Rev. D 101 , 122001 (2020), arXiv:1908.06050 [grqc].\n- [58] H.-Y. Chen, M. Fishbach, and D. E. Holz, A two per cent Hubble constant measurement from standard sirens within five years, Nature 562 , 545 (2018), arXiv:1712.06531 [astro-ph.CO].\n- [59] S. Mastrogiovanni, D. Laghi, R. Gray, G. C. Santoro, A. Ghosh, C. Karathanasis, K. Leyde, D. A. Steer, S. Perriès, and G. Pierra, Joint population and cosmological properties inference with gravitational waves standard sirens and galaxy surveys, Phys. Rev. D 108 , 042002 (2023), arXiv:2305.10488 [astro-ph.CO].\n- [60] N. Borghi, M. Mancarella, M. Moresco, M. Tagliazucchi, F. Iacovelli, A. Cimatti, and M. Maggiore, Cosmology and Astrophysics with Standard Sirens and Galaxy Catalogs in View of Future Gravitational Wave Observations, Astrophys. J. 964 , 191 (2024), arXiv:2312.05302 [astro-ph.CO].\n- [61] R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. X. Adhikari, V. B. Adya, C. Affeldt, D. Agarwal, M. Agathos, K. Agatsuma, N. Aggarwal, O. D. Aguiar, L. Aiello, A. Ain, P. Ajith, S. Akcay, T. Akutsu, S. Albanesi, A. Allocca, P. A.  \n- Altin, A. Amato, C. Anand, S. Anand, A. Ananyeva, S. B. Anderson, and Others (LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration), Gwtc-3: Compact binary coalescences observed by ligo and virgo during the second part of the third observing run, Phys. Rev. X 13 , 041039 (2023).\n- [62] A. M. Farah, B. Edelman, M. Zevin, M. Fishbach, J. M. Ezquiaga, B. Farr, and D. E. Holz, Things that might go bump in the night: Assessing structure in the binary black hole mass spectrum, The Astrophysical Journal 955 , 107 (2023).\n- [63] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, One channel to rule them all? constraining the origins of binary black holes using multiple formation pathways, The Astrophysical Journal 910 , 152 (2021).\n- [64] DESI Collaboration, The desi experiment part i: Science,targeting, and survey design (2016), arXiv:1611.00036 [astro-ph.IM].\n- [65] D. J. Schlegel, S. Ferraro, G. Aldering, C. Baltay, S. BenZvi, R. Besuner, G. A. Blanc, A. S. Bolton, A. Bonaca, D. Brooks, E. Buckley-Geer, Z. Cai, J. DeRose, A. Dey, P. Doel, A. Drlica-Wagner, X. Fan, G. Gutierrez, D. Green, J. Guy, D. Huterer, L. Infante, P. Jelinsky, D. Karagiannis, S. M. Kent, A. G. Kim, J.-P. Kneib, J. A. Kollmeier, A. Kremin, O. Lahav, M. Landriau, D. Lang, A. Leauthaud, M. E. Levi, E. V. Linder, C. Magneville, P. Martini, P. McDonald, C. J. Miller, A. D. Myers, J. A. Newman, P. E. Nugent, N. Palanque-Delabrouille, N. Padmanabhan, A. Palmese, C. Poppett, J. X. Prochaska, A. Raichoor, S. Ramirez, N. Sailer, E. Schaan, M. Schubnell, U. Seljak, H.-J. Seo, J. Silber, J. D. Simon, Z. Slepian, M. Soares-Santos, G. Tarle, M. Valluri, N. J. Weaverdyck, R. H. Wechsler, M. White, C. Yeche, and R. Zhou, A spectroscopic road map for cosmic frontier: Desi, desi-ii, stage-5 (2022), arXiv:2209.03585 [astro-ph.CO].\n- [66] K. E. S. Ford and B. McKernan, Ligo tells us liners are not optically thick riafs, Monthly Notices of the Royal Astronomical Society: Letters 490 , L42-L46 (2019)."}