© The Authors 2025

1. Introduction

The first direct detection of gravitational waves (GWs) with GW150914 (Abbott et al. 2016a) marked the beginning of a new era in astronomy and cosmology. GWs from compact binary coalescences (CBC) can be used as “standard sirens” as they provide a direct measurement of the luminosity distance to the source, without the need for any intermediate calibrator (Schutz 1986; Holz & Hughes 2005; Moresco et al. 2022). By combining standard sirens with information on the redshift of the source, it becomes possible to study the expansion history of the Universe through the electromagnetic luminosity distance-redshift relation:

$$\begin{array}{c}\hfill d_L^{\mathrm{em}}(z;{\mathit{\lambda }}_c)=(1+z){\displaystyle {\int }_0^z}\frac{c\mathrm{d}z}{H(z;{\mathit{\lambda }}_c)},\end{array}$$(1)

where we assumed a flat geometry of the Universe, and λ_c represents cosmological parameters such as the Hubble constant, H₀, and the matter energy density, Ω_m. Standard sirens are, and will increasingly become, a powerful tool for addressing tensions in cosmological parameters, such as the Hubble tension, that have emerged in the recent era of precision cosmology (Verde et al. 2019; Moresco et al. 2022). Interestingly, standard sirens also allow us to constrain possible deviations in the propagation of gravitational and electromagnetic signals, providing an additional test of general relativity (GR). The multimessenger observation of the binary neutron star coalescence GW170817 (Abbott et al. 2017a,b) placed a stringent limit on the difference between the propagation speed of gravitational and electromagnetic waves, |c_gw − c|/c < 𝒪(10⁻¹⁵), ruling out many modified-gravity (MG) models. All MG models consistent with this constraint introduce a friction term in the propagation of tensor modes at cosmological scales, which is absent in GR (Belgacem et al. 2019a). This friction term modifies the effective distance to which standard sirens are sensitive, making them an effective tool to probe MG models (Lombriser & Taylor 2016; Saltas et al. 2014; Nishizawa 2018; Arai & Nishizawa 2018; Amendola et al. 2018; Belgacem et al. 2018a,b). A common way to parameterize the friction term involves two parameters λ_mg = (Ξ₀, n). In this parameterization, the effective luminosity distance measured by standard sirens is given by (Belgacem et al. 2018a)

$$\begin{array}{cc}\hfill d_L^{\mathrm{gw}}(z;{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}})& =d_L^{\mathrm{em}}(z;{\mathit{\lambda }}_c)\mathrm{\Xi }(z;{\mathit{\lambda }}_{\mathrm{mg}})\hfill \\ \hfill & =d_L^{\mathrm{em}}(z;{\mathit{\lambda }}_c)[{\mathrm{\Xi }}_0+\frac{1-{\mathrm{\Xi }}_0}{{(1+z)}^n}]\hfill \end{array}$$(2)

where Ξ(z; λ_mg) quantifies the deviation from the standard electromagnetic distance Eq. (1) at each redshift. This parameterization can represent several scalar-tensor theories of the Horndeski class, including the Brans-Dicke, f(R), covariant-Galileon, and minimal-acceleration models (see Table 1 in Belgacem et al. 2018b for a summary). It can also be connected with nonlocal gravity theories (Belgacem et al. 2019b, 2020).

One of the main challenges to exploiting the standard siren method for cosmological purposes is the perfect degeneracy between the binary redshift and the chirp mass. This degeneracy, due to the scale invariance of GR (Mastrogiovanni & Steer 2022; Mastrogiovanni et al. 2024a), can be broken if a physical effect imprints a known scale on the GW signal or if external redshift information is provided. An example of a standard siren technique that determines the redshift by leveraging physical scales in the GW signal is the “spectral-siren” method, which exploits features in the source-frame mass distribution of CBC populations to statistically infer the redshift (Chernoff & Finn 1993; Taylor et al. 2012; Farr et al. 2019; Mastrogiovanni et al. 2021; Mukherjee 2022; Mancarella et al. 2022; Karathanasis et al. 2023; Chen et al. 2024a). Examples of such features include those observed in the binary black hole (BBH) mass distribution inferred from LIGO-Virgo-KAGRA (LVK) data up to the observing run O3, which reveals an overdensity at approximately 35 M_⊙ and a steep decrease after 80 M_⊙ (Abbott et al. 2023a). There are two main ways to include external redshift information in the inference process. The most intuitive is to use, if identified, the redshift of the electromagnetic counterpart when the standard siren is “bright” (Holz & Hughes 2005; Nissanke et al. 2010); this requires the presence of an electromagnetic phenomenon associated with the CBC, such as a kilonova with the merger of a binary neutron star. The other method, also known as the “galaxy catalog method”, requires the GW luminosity-distance probability distribution to be statistically combined with a redshift prior distribution constructed from a catalog of potential hosts, usually a galaxy catalog (Schutz 1986; Del Pozzo 2012; Fishbach et al. 2019; Gray et al. 2020; Gair et al. 2023). All three approaches rely on assumptions about the underlying population, whether it is the host galaxy for bright sirens, the redshift prior in the galaxy-catalog method, or the mass distribution in the spectral-siren case. Various pipelines have been developed to implement the spectral-siren method, namely POPMODELS (Wysocki & O’Shaughnessy 2017), GWPopulation (Talbot et al. 2019), MGCosmoPop (Mancarella et al. 2022), ICAROGW 1.0 (Mastrogiovanni et al. 2021), SODAPOP (Landry 2021), GWInferno (Edelman et al. 2023), and in Chen et al. (2024a). More recently, a differentiable pipeline that samples the full hierarchical population posterior was released (Mancarella & Gerosa 2025). For the galaxy-catalog method, tools such as DarkSirensStat (Finke et al. 2021), GWCOSMO (Gray et al. 2023), and – specifically for LISA – cosmolisa (Laghi et al. 2021) are available. In recent years, new codes have been developed that unify the different methods in a unique Bayesian framework, jointly inferring cosmological and MG hyperparameters with population ones to provide robust constraints: ICAROGW 2.0 (Mastrogiovanni et al. 2023, 2024b), GWCOSMO (Gray et al. 2023), and CHIMERA (Borghi et al. 2024) (referred to as CHIMERA 1.0 from now on).

Constraints on the Hubble constant H₀ from standard sirens have been obtained using the publicly available GWTC-3 data (Abbott et al. 2019, 2021a, 2023b). These include measurements from the bright siren GW170817 (Abbott et al. 2017c; Palmese et al. 2024), as well as dark sirens in GWTC-3 combined with the GLADE+ (Dálya et al. 2022) galaxy catalog (Fishbach et al. 2019; Abbott et al. 2021b, 2023c; Finke et al. 2021; Mastrogiovanni et al. 2023; Gray et al. 2023), the DES survey (Soares-Santos et al. 2019; Palmese et al. 2020), the DELVE survey Alfradique et al. (2024), or DESI (Ballard et al. 2023). More recent work has extended these analyses to dark sirens from the O4 run (Bom et al. 2024). Similarly, constraints on Ξ₀ have been derived from GWTC-3 data. For example, Finke et al. (2021) derived ${\mathrm{\Xi }}_0=2.1_{-2.1}^{+3.2}$ using dark sirens and GLADE+, while Mancarella et al. (2022) constrained Ξ₀ with 58% uncertainty using the spectral siren method. More recently, joint constraints on H₀ and Ξ₀ were obtained by combining the spectral siren and galaxy catalog methods, yielding am uncertainty of 73% on Ξ₀ and 58% on H₀ using 42 BBH events of GWTC-3 and the GLADE+ galaxy catalog (Chen et al. 2024b). Preliminary forecasts on future constraints for MG and cosmological parameters have been explored using spectral sirens. For instance, Leyde et al. (2022) predicts a 20 − 30% measurement of Ξ₀ using about 400 detections at the design LVK O5 sensitivity in a GR scenario. Forecasts using bright sirens (e.g., Niu et al. 2021; Chen et al. 2024c; Colangeli et al. 2025) or alternative approaches, such as the GW galaxies cross-correlation (e.g., Mukherjee et al. 2021; Afroz & Mukherjee 2024a,b), have also been investigated.

In general, standard siren codes are computationally limited by the number of GW events they can process. To forecast how next-generation interferometers, such as the Einstein Telescope (Branchesi et al. 2023; Abac et al. 2025), Cosmic Explorer (Reitze et al. 2019), and LISA (Colpi et al. 2024), will constrain cosmology and modified GW propagation, it is necessary to improve existing pipelines, as these experiments will detect up to 10⁵ BBHs per year.

This paper has two primary aims. First, we present CHIMERA 2.0, an enhanced version of CHIMERA 1.0 that can handle up to tens of thousands of GW events, a critical step toward next-generation detector data. This is achieved by introducing three different kernel-density-estimate (KDE) algorithms, which form the backbone of the pipeline. These KDEs provide better flexibility and fully leverage GPU acceleration for optimal performance. In a future study, we will test and validate this code against other existing pipelines through a blinded-mock-data challenge. This work will assess how the computational cost of current codes scales with the number of events and identify potential systematic effects in the implemented algorithms. Second, we used this upgraded code, that includes MG models, to forecast joint cosmological and MG constraints for the first time for the future O5 observing run of the LVK collaboration using the combination of spectral-siren and galaxy-catalog methods.

This paper is organized as follows. Section 2 outlines the statistical framework of combined standard siren methods and the enhanced numerical implementation of CHIMERA 2.0, comparing it with CHIMERA 1.0 in terms of results and computational efficiency. Section 3 describes the mock catalogs used to forecast O5-like constraints on cosmology and modified GW propagation. We studied three different scenarios: one in which there is no MG, one with a Ξ₀ greater than 1, and one with Ξ₀ < 1. Finally, Section 4 presents the results.

2. Methods

2.1. Statistical framework

Standard sirens as cosmological probes provide constraints on the hyperparameters (Λ) describing the properties of the CBC population and the underlying cosmological model from a set of observations drawn from that population. Observations are incomplete due to selection biases in GW interferometers and noisy because the single event parameters in detector-frame {θ_i^d}_i = 1^N_obs (e.g., binary redshifted masses, spins, luminosity distance, and localization area) are inferred from a set of observations {d_i}_i = 1^N_obs including measurement uncertainties. The correct Bayesian framework that accounts for selection effects and noisy observations is a hierarchical inhomogeneous Poisson process, described by the hyper-likelihood (Loredo 2004; Mandel et al. 2019; Vitale et al. 2020):

$$\begin{array}{cc}\hfill \mathcal{L}({\{{\mathit{d}}_i\}}_{i=1}^{N_{\mathrm{obs}}}\mid \mathbf{\Lambda })& \propto \frac{1}{{[\xi (\mathbf{\Lambda })]}^{N_{\mathrm{obs}}}}{\displaystyle \prod _{i=1}^{N_{\mathrm{obs}}}\int \mathrm{d}}{\mathit{\theta }}_i^{\mathrm{d}}\frac{p_{\mathrm{gw}}({\mathit{\theta }}_i^{\mathrm{d}}\mid {\mathit{d}}_i)}{\pi ({\mathit{\theta }}_i^{\mathrm{d}})}{p}_{\mathrm{pop}}({\mathit{\theta }}_i^{\mathrm{d}}\mid \mathbf{\Lambda }),\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill \xi (\mathbf{\Lambda })& ={\displaystyle \int \mathrm{d}}{\mathit{\theta }}^{\mathrm{d}}{P}_{\det }({\mathit{\theta }}^{\mathrm{d}})p_{\mathrm{pop}}({\mathit{\theta }}^{\mathrm{d}}\mid \mathbf{\Lambda }).\hfill \end{array}$$(4)

ξ(Λ), representing the fraction of population that can be detected, corrects the Malmquist bias due to selection effects and depends on the probability, P_det(θ^d), of detecting a source with detector-frame parameters θ^d:

$$\begin{array}{c}\hfill P_{\det }\left({\mathit{\theta }}^{\mathrm{d}}\right)={\displaystyle {\int }_{\mathit{d}\in \text{detectable}}\mathrm{d}\mathit{d}}\frac{p_{\mathrm{gw}}({\mathit{\theta }}^{\mathrm{d}}\mid \mathit{d})}{\pi ({\mathit{\theta }}^{\mathrm{d}})}·\end{array}$$(5)

The term p_gw is the posterior distribution for the detector-frame parameters given the data, and it accounts for the measurement uncertainties. The prior distribution, π, which appears in the denominator of the integrand, converts this posterior to its corresponding likelihood. The population prior, p_pop, gives the probability of drawing an event with parameters θ^d from a population described by hyperparameters, Λ.

The term p_pop is usually modeled in terms of source frame parameters θ. Here, we neglected spins and focused only on binary masses, m_1, 2, redshift, z, and sky position, $\widehat{\mathrm{\Omega }}$. The cosmological and MG hyperparameter (λ_c, λ_mg) map between detector- and source-frame parameters. Specifically, they convert the measured luminosity distance into the corresponding redshift via Eq. (2). The redshift is then used to transform the binary masses in the source frame as m_1, 2^d = (1 + z) m_1, 2. By assuming that the mass distribution does not evolve with redshift, the source-frame population prior can be factorized as

$$\begin{array}{c}\hfill p_{\mathrm{pop}}(\mathit{\theta }\mid \mathbf{\Lambda })=p(m_1,m_2\mid {\mathit{\lambda }}_m)p_{\mathrm{cbc}}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{\lambda }}_c,{\mathit{\lambda }}_r),\end{array}$$(6)

where the additional hyperparameters λ_m and λ_r describe the mass and merger-rate evolution distributions, respectively. The set of all different hyperparameters is denoted by Λ = {λ_c, λ_mg, λ_m, λ_r}. The first term in Eq. (6) describes the mass distribution of CBCs, while p_cbc is the probability of having a CBC at redshift, z, and sky position $\widehat{\mathrm{\Omega }}$. The latter probability can be written as the probability of having a host galaxy at $(z,\widehat{\mathrm{\Omega }})$, denoted by $p_{\mathrm{gal}}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{\lambda }}_c)$, times a function describing the merger-rate redshift evolution, denoted by ψ(z; λ_r). The distribution p_gal is expressed as (Chen et al. 2018; Finke et al. 2021):

$$\begin{array}{c}\hfill p_{\mathrm{gal}}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{\lambda }}_c)=f_R{p}_{\mathrm{cat}}(z,\widehat{\mathrm{\Omega }})+(1-f_R)p_{\mathrm{miss}}(z,\widehat{\mathrm{\Omega }}),\end{array}$$(7)

where

$$\begin{array}{cc}\hfill p_{\mathrm{cat}}(z,\widehat{\mathrm{\Omega }})& \propto {\displaystyle \sum _g}{w}_g\delta (\mathrm{\Omega }-{\mathrm{\Omega }}_g)\frac{\mathcal{N}(z;{\stackrel{\sim }{z}}_g,{\stackrel{\sim }{\sigma }}_{z,g})\frac{\mathrm{d}{V}_c}{\mathrm{d}z}(z;{\mathit{\lambda }}_c)}{{\displaystyle \int \mathrm{d}}z\mathcal{N}(z;{\stackrel{\sim }{z}}_g,{\stackrel{\sim }{\sigma }}_{z,g})\frac{\mathrm{d}{V}_c}{\mathrm{d}z}(z;{\mathit{\lambda }}_c)},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill p_{\mathrm{miss}}(z,\widehat{\mathrm{\Omega }})& =\frac{1-P_{\mathrm{compl}}(z,\widehat{\mathrm{\Omega }})}{1-f_R}\frac{1}{V_c(z_{\max };{\mathit{\lambda }}_c)}\frac{\mathrm{d}{V}_c}{\mathrm{d}z}(z;{\mathit{\lambda }}_c),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill f_R& =\frac{1}{V_c(z_{\max };{\mathit{\lambda }}_c)}{\displaystyle \int \mathrm{d}}{V}_c{P}_{\mathrm{compl}}(z,\widehat{\mathrm{\Omega }}).\hfill \end{array}$$(10)

In this framework, p_cat is the probability distribution derived from the galaxy catalog and is modeled as a sum of Gaussian distributions centered on the measured galaxy redshifts (${\stackrel{\sim }{z}}_g$), with standard deviations equal to the measurement uncertainties (${\stackrel{\sim }{\sigma }}_g)$. Each Gaussian is multiplied by a prior distribution, assumed to be uniform in comoving volume in the absence of additional information (Gair et al. 2023). The galaxy positions on the sky are assumed to be known without uncertainty. The term p_miss accounts for the missing galaxies, and depends on P_compl, the probability of missing a galaxy at $(z,\widehat{\mathrm{\Omega }})$, as well as on assumptions on how missing galaxies are distributed. In Eq. (9), they are assumed to be homogeneously distributed. However, more accurate prescriptions that account for galaxy clustering properties have recently been proposed (Finke et al. 2021; Dalang & Baker 2024; Leyde et al. 2024; Dalang et al. 2024). Finally, f_R (10) is the galaxy-catalog completeness fraction.

To use the source-frame population prior (6) in Eq. (3) it is necessary to take into account the Jacobian of the transformation θ^d → θ(θ^d; λ_c, λ_mg):

$$\begin{array}{c}\hfill p_{\mathrm{pop}}({\mathit{\theta }}^{\mathrm{d}}\mid \mathbf{\Lambda })\propto p_{\mathrm{pop}}(\mathit{\theta }\mid \mathbf{\Lambda })\times {\left|\frac{\mathrm{d}{\mathit{\theta }}^{\mathrm{d}}}{\mathrm{d}\mathit{\theta }}\right|}^{-1}=\frac{p_{\mathrm{pop}}(\mathit{\theta }\mid \mathbf{\Lambda })}{{(1+z)}^3\left|\frac{\partial d_L}{\partial z}(z;{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}})\right|},\end{array}$$(11)

where one (1 + z) factor expresses the conversion of the merger rate from the source to the detector-frame, and the other two come from the redshifting of binary masses, and the term |∂d_L/∂z| is due to the luminosity distance-redshift conversion.

2.2. chimera1.0

In this section, we describe the numerical implementation of CHIMERA 1.0 and its main computational bottlenecks.

The term related to the selection bias ξ(Λ) (4) is estimated using a set of N_inj simulated events (injections) with parameters, {θ_j^d}_j = 1^N_inj, that span the detectable parameter space. The integral in Eq. (4) is approximated by the following Monte Carlo summation over the set of injections (Talbot et al. 2019; Thrane & Talbot 2019; Essick & Farr 2022):

$$\begin{array}{c}\hfill \xi (\mathbf{\Lambda })\approx \frac{1}{N_{\mathrm{inj}}}{\displaystyle \sum _{j=1}^{N_{\det }}}\frac{p_{\mathrm{pop}}({\mathit{\theta }}_j\mid \mathbf{\Lambda })}{{(1+z_j)}^3\left|\frac{\partial d_L}{\partial z}(z_j;{\mathit{\lambda }}_{\mathit{c}},{\mathit{\lambda }}_{\mathrm{mg}})\right|}\equiv \frac{1}{N_{\mathrm{inj}}}{\displaystyle \sum _{j=1}^{N_{\det }}}{s}_j,\end{array}$$(12)

where we inserted Eqs. (5) and (11) into Eq. (4). We checked the numerical stability of the previous finite Monte Carlo summation by requiring that the “effective” number of injections, defined as in Farr (2019)

$$\begin{array}{c}\hfill N_{\mathrm{eff}}^{\mathrm{inj}}={\left[{\displaystyle \sum _{j=0}^{N_{\det }}}{s}_j\right]}^2\times {[{\displaystyle \sum _{j=0}^{N_{\det }}}{s}_j^2-\frac{1}{N_{\mathrm{inj}}}{\left({\displaystyle \sum _{j=0}^{N_{\det }}}{s}_j\right)}^2]}^{-1},\end{array}$$(13)

is larger than 5N_det.

The N_obs integrals appearing in the numerator of Eq. (3), which we denote as I_i(d_i ∣ Λ), are calculated in the source frame over the three-dimensional volume defined by the GW event localization area, $\delta {\widehat{\mathrm{\Omega }}}_i$, and redshift interval, δz_i(λ_c, λ_mg):

$$\begin{array}{cc}& I_i({\mathit{d}}_i\mid \mathbf{\Lambda })={\int }_{\delta {\mathrm{\Omega }}_i\times \delta z_i({\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}})}{\mathrm{d}}^2\widehat{\mathrm{\Omega }}\mathrm{d}z\hfill \\ \hfill & \times {\mathcal{K}}_{\mathrm{gw},i}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{d}}_i,{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}},{\mathit{\lambda }}_m)\frac{p_{\mathrm{cbc}}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{\lambda }}_c,{\mathit{\lambda }}_r)}{{(1+z)}^3\left|\frac{\partial d_L}{\partial z}(z;{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}})\right|}.\hfill \end{array}$$(14)

Here, we used the transformation dθ_i^d p_gw(θ_i^d ∣ d_i) = dθ_i p_gw(θ_i ∣ d_i; λ_c, λ_mg). The integrand in Eq. (14) includes terms from the population prior (6) and the Jacobian (11) that depend only on the redshift and/or sky position. These terms are multiplied by the GW event kernel, 𝒦_gw, i, defined as the GW posterior weighted by the mass distribution and detector-frame parameter priors, marginalized over m_1, 2 and evaluated on the integration volume $(z,\widehat{\mathrm{\Omega }})$:

$$\begin{array}{cc}& {\mathcal{K}}_{\mathrm{gw},i}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{d}}_i,{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}},{\mathit{\lambda }}_m)={\displaystyle \int \mathrm{d}}{m}_{1,i}\mathrm{d}{m}_{2,i}\hfill \\ \hfill & \times {p_{\mathrm{gw}}(m_{1,i},m_{2,i},z_i,{\widehat{\mathrm{\Omega }}}_i\mid {\mathit{d}}_i;{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}})\frac{p(m_{1,i},m_{2,i}\mid {\mathit{\lambda }}_m)}{\pi (m_{1,i}^{\mathrm{d}},m_{2,i}^{\mathrm{d}})\pi (d_{L,i})}|}_{(z,\widehat{\mathrm{\Omega }})}.\hfill \end{array}$$(15)

Here, we discuss the “3D kernel”. The GW kernel (15) is approximated using a weighted KDE built using N_s posterior estimate (PE) samples drawn from p_gw. This algorithm employs a Gaussian kernel and a 3D training dataset:

$$\begin{array}{cc}& {\mathcal{K}}_{\mathrm{gw},i}(z,\widehat{\mathrm{\Omega }}\mid {\mathit{d}}_i,{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}},{\mathit{\lambda }}_m)\hfill \\ \hfill & \approx {\text{KDE}\left[\{(z_i^j,{\widehat{\mathrm{\Omega }}}_i^j){|w_i^j=\frac{p(m_{1,i}^j,m_{2,i}^j\mid {\mathit{\lambda }}_m)}{\pi (m_{1,i}^{\mathrm{d},j},m_{2,i}^{\mathrm{d},j})\pi (d_{L,i}^j)}\}}_{j=1}^{N_{\mathrm{s}}}\right]|}_{(z,\widehat{\mathrm{\Omega }})}.\hfill \end{array}$$(16)

Here, the index j refers to the PE samples of the i-th event. This KDE is evaluated on a 3D grid that approximates the integration domain of Eq. (14); it is constructed as follows:

1. Divide the 2D localization area, $\delta {\widehat{\mathrm{\Omega }}}_i$, into N_pixⁱ equal-area pixels (see top panel of Fig. 1), as first proposed in Gray et al. (2022).

   Fig. 1. Top: Pixelization of 90% localization area of a mock GW event, with PE samples (marked with a cross) colored according to the pixel they fall into. Bottom: Different GW kernels implemented in CHIMERA 2.0. Each line represents the GW kernel in a pixel of the above mock GW event. Kernels are computed at single hyperparameter space points.

2. Discretize the redshift interval, δz_i(λ_c, λ_mg), into N_z equally spaced points, ensuring that the grid covers all possible redshifts given the prior on cosmological and MG hyperparameters. This allows us to significantly reduce the computational cost by computing the p_cat term only once before the inference.

3. Repeat the redshift grid N_pix times and duplicate the pixel center coordinates N_z times to match the redshifts within each pixel.

In CHIMERA 1.0, the KDE evaluation is the main computational bottleneck, accounting for approximately 80% of the total time required for a complete population fit in a scenario involving 100 GW events – with 5000 PE samples per event – and 2 × 10⁷ injections. The cumulative computational time for the KDE evaluation scales as follows:

$$\begin{array}{c}\hfill t_{\mathrm{KDE}}\sim \mathcal{O}(N_{\mathrm{obs}}\times N_{\mathrm{z}}\times N_{\mathrm{pix}}\times N_{\mathrm{s}}),\end{array}$$(17)

where N_z, N_pix are the resolutions of redshift and localization area integration volumes, respectively. In the above scenario, the population fit required 10⁴ CPU hours using the emcee sampler (Foreman-Mackey et al. 2013) with 50 walkers parallelized across 25 Intel CPUs (2.0 GHz, 1 core per CPU) on a HPC facility. Due to the linear dependence on the number of GW events, the KDE bottleneck imposes a significant computational challenge for future next-generation interferometers that will detect hundreds of thousands of GW events. For example, CHIMERA 1.0 analyzes around 1000 events in more than a month of CPU time – a feasible but increasingly impractical burden for even larger catalogs. This computational time cannot be reduced by simply allocating more CPUs due to inherent limitations in most of the sampling algorithms. Instead, more efficient algorithms and hardware accelerators as the GPU were preferred in order to reduce the computational burden of the likelihood evaluation.

2.3. Enhanced numerical implementation

To overcome the computational limitations of CHIMERA 1.0, we introduced a new release featuring a redesigned pipeline architecture and advanced KDE algorithms.

The new code release, CHIMERA 2.0¹, is fully implemented in the JAX framework (Bradbury et al. 2018). JAX is a high-performance Python library for numerical computing and machine learning that combines NumPy-like syntax with just-in-time compilation, automatic differentiation, and GPU or TPU acceleration. This enables CHIMERA 2.0 to leverage gradient-based Markov chain Monte Carlo (MCMC) algorithms, such as Hamiltonian MCMC algorithms. Population and cosmological models are constructed using equinox (Kidger & Garcia 2021), which ensures data structures are compatible as JAX pytrees and supports hyperparameter vectorization. Additionally, CHIMERA 2.0 now includes extended cosmological models, such as curved and evolving dark-energy universes.

CHIMERA 2.0 includes three different KDE algorithms, which have been tested and validated against each other. The first is the “3D” kernel implemented in the first release. However, in CHIMERA 2.0, it is evaluated only on the portion of the redshift grid containing redshift samples. Since the redshift grids must be pre-computed considering the cosmological priors, the redshift samples are contained in a smaller fraction – up to 1/3 – of the redshift grids at each MCMC step. In the remaining segments of the grids, the kernel is set to zero.

In the “many-1D” kernel approach, the pixelization and the pre-computation of the redshift grids and p_cat remain the same as in the previous case. However, this algorithm approximates the GW kernel in each pixel of δΩ_i, indexed with k, using 1D weighted KDEs for each pixel, for a total of N_pixⁱ KDEs. Each KDE is trained using the N_kⁱ PE samples that are comprised in the corresponding pixel, as illustrated in the top panel of Fig. 1. In other words, this algorithm computes the 1D redshift distributions for each pixel by marginalizing over the sky localization distribution. The GW kernel within each pixel is then obtained by multiplying the 1D KDE by a pixel-dependent normalization factor given by the 2D KDE of ${\{\widehat{\mathrm{\Omega }}\}}_{j=1}^{N_{\mathrm{s}}}$ evaluated at the center of the pixel:

$$\begin{array}{cc}& {\mathcal{K}}_{\mathrm{gw},i}(z,{\widehat{\mathrm{\Omega }}}_k\mid {\mathit{d}}_i,{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}},{\mathit{\lambda }}_m)\hfill \\ \hfill & \approx {\text{KDE}\left[\{z_i^j{|w_i^j=\frac{p(m_{1,i}^j,m_{2,i}^j\mid {\mathit{\lambda }}_m)}{\pi (m_{1,i}^{\mathrm{d},j},m_{2,i}^{\mathrm{d},j})\pi (d_{L,i}^j)}\}}_{j=1}^{N_k^i}\right]|}_z\hfill \\ \hfill & \times {\text{KDE}\left[{\left\{{\widehat{\mathrm{\Omega }}}_i^j\right\}}_{j=1}^{N_{\mathrm{s}}}\right]|}_{{\widehat{\mathrm{\Omega }}}_k}.\hfill \end{array}$$(18)

The normalization factors can be pre-computed once, as they are independent of hyperparameters. Also in this case, at each MCMC step, the KDEs are evaluated only on the portions of the pre-computed redshift grids containing redshift samples. Two kernel options are available for the 1D KDEs: Gaussian and Epanechnikov, with the latter being slightly faster. To ensure a consistent dataset dimension across pixels and events, the weighted PE samples in each pixel are binned. This enables vectorized computation over both pixels and events, eliminating the costly Python loops used in the previous algorithm.

The “single-1D” kernel algorithm is very similar to the “many-1D” approach, but only requires a single KDE computation per event rather than N_pixⁱ KDEs. Instead of splitting the samples by pixel, all N_s redshift samples are used to train a weighted 1D KDE, which represents the redshift distribution of the GW event weighted by the mass distribution. The GW kernel in each pixel is then estimated by multiplying this global weighted redshift distribution by the same pixel-dependent normalization factors:

$$\begin{array}{cc}& {\mathcal{K}}_{\mathrm{gw},i}(z,{\widehat{\mathrm{\Omega }}}_k\mid {\mathit{d}}_i,{\mathit{\lambda }}_c,{\mathit{\lambda }}_{\mathrm{mg}},{\mathit{\lambda }}_m)\hfill \\ \hfill & \approx {\text{KDE}\left[\{z_i^j{|w_i^j=\frac{p(m_{1,i}^j,m_{2,i}^j\mid {\mathit{\lambda }}_m)}{\pi (m_{1,i}^{\mathrm{d},j},m_{2,i}^{\mathrm{d},j})\pi (d_{L,i}^j)}\}}_{j=1}^{N_{\mathrm{s}}}\right]|}_z\hfill \\ \hfill & \times {\text{KDE}\left[{\left\{{\widehat{\mathrm{\Omega }}}_i^j\right\}}_{j=1}^{N_{\mathrm{s}}}\right]|}_{{\widehat{\mathrm{\Omega }}}_k}.\hfill \end{array}$$(19)

This algorithm assumes independence between the sky position and redshift. Thus, this approximation is suitable when these two variables are weakly correlated. As in the “many-1D” case, this KDE supports Gaussian or Epanechnikov kernels, its evaluation is restricted to the relevant portion of the pre-computed redshift grid, and the dataset can be binned. Also, this algorithm is fully vectorized over the events.

2.4. Scaling performances

The bottom panels of Fig. 1 compare the three different GW kernels considering a mock GW event with the localization area divided into 8 pixels. The “3D” kernel distributions show fluctuations in less populated pixels. In contrast, the “single-1D” kernel smooths out these substructures, while the “many-1D” one provides only partial smoothing. In Appendix A we compare and validate the “single-1D” and “many-1D” against the implementation of CHIMERA 1.0. In particular, Fig. A.1 highlights that the three kernels produce results that are consistent with each other, with no significant differences observed.

The new kernels yield the likelihood evaluation times shown in the top panel of Fig. 2. Since the “3D” kernel is only calculated on a limited portion of the pre-computed redshift grid, the dependence of t_KDE on N_z is reduced, achieving a 2 − 3 speed improvement over CHIMERA 1.0. However, this algorithm still relies on a computationally expensive loop over the events, resulting in a linear increase of t_KDE with N_obs. The vectorized “many-1D” and “single-1D” algorithms mitigate this linear scaling, particularly on a GPU. The computational time for the “many-1D” kernel is further reduced compared to the “3D” approach due to the dataset’s lower dimensionality, resulting in a 10 − 20 (100 − 200) speed improvement over CHIMERA 1.0 on the CPU (GPU), depending on the number of events. The “single-1D” case, in turn, does not depend on the number of pixels and benefits from binning to reduce the dependence of t_KDE on N_s, making it the fastest algorithm among the three. This algorithm is 20 − 100 (100 − 1000) times faster than the first release on the CPU (GPU), depending on the number of events.

Fig. 2. Top: Likelihood evaluation times for catalogs with different numbers of events (each with 5000 PE samples) at a single hyperparameter space point, compared across different kernel algorithms. Solid curves show times on four Intel CPUs (2.0 GHz clock frequency), while dashed curves represent times on a single Nvidia A100 GPU. Bottom: Run times for a full affine-invariant MCMC fit, using the emcee sampler with 50 walkers, on 25 CPUs (left panel) or one GPU (right panel).

The “many-1D” kernel reduces the time spent on KDE analysis from 80% (as in CHIMERA 1.0) to 35% for a fit of a catalog of 100 GW events with 5000 PE samples and 2 × 10⁷ injections. The “single-1D” algorithm further cuts this fraction to 15%. With the “many-1D” algorithm, the GW kernel (15) and bias term (4) evaluations take a similar time, but for larger datasets, the kernel becomes again the main bottleneck. In contrast, in the “single-1D” case, the total time is dominated by the bias term, which remains the limiting factor even with more events – since the required injections for numerical accuracy is expected to scale as ∝N_obs² (Talbot & Golomb 2023).

The bottom panels of Fig. 2 show the total run time of a full affine-invariant MCMC fit using the emcee sampler with 50 walkers. In the left panel, run times are shown for 25 CPUs, the maximum number that can be used to parallelize walker moves. In the right panel, run times on a single GPU are presented, where walkers evolve sequentially. Poorly vectorized algorithms, such as CHIMERA 1.0 and the “3D” kernel, are less efficient on GPUs than on multiple CPUs, despite faster single-point evaluation on GPUs. In contrast, the “many-1D” and “single-1D” kernels perform better on GPUs and exhibit a weaker linear dependence on the number of events. These KDE algorithms are a key advancement, preparing CHIMERA 2.0 for next-generation detectors that will detect 1000 times more events than current observatories. Finally, we stress that for less parallelizable methods, such as standard MCMC, parallel tempering MCMC, Hamiltonian MCMC, or nested-sampling methods – where fewer chains are evolved in parallel compared to affine-invariant MCMC options – the advantage of a single GPU over multiple CPUs becomes even more significant.

Although the “single-1D” kernel is the most efficient and produces results consistent with the other algorithms, CHIMERA 2.0 retains all three KDE methods. This ensures flexibility, allowing the code to handle also real GW events with complex posterior distributions, where the assumption behind the “single-1D” kernel may not hold.

3. Data

In the following, we outline the process used to generate the three mock GW catalogs for the MCMC analyses, following a method similar to that of Borghi et al. (2024).

3.1. GW populations

The mock GW catalogs were generated starting from a mock galaxy catalog and populating galaxies with CBC events drawn from a fiducial population model. These are detailed in the following paragraphs.

Similarly to Borghi et al. (2024), the galaxy catalog considered in this work is a complete subsample from the MICE Grand Challenge light-cone simulation (v2) (Carretero et al. 2015); Fosalba:2013wxa; Fosalba:2013mra; Hoffmann:2014ida. The MICEv2 catalog covers one octant of the sky and is designed to mimic a complete DES-like survey up to an observed magnitude of i < 24 at redshift z < 1.4. The galaxy catalog was obtained by considering only galaxies with stellar masses log M_*/M_⊙ > 10.5 and with a uniform comoving volume redshift distribution up to z < 1.3, resulting in 1.6 million galaxies. This cut is aligned with the assumption that the binary merger rate follows the stellar mass distribution. This assumption is widely adopted in the current literature through absolute magnitude cuts and luminosity weighting (Fishbach et al. 2019; Finke et al. 2021; Gray et al. 2022; Abbott et al. 2023c; Muttoni et al. 2023; Alfradique et al. 2025). In the left panel of Fig. 3, we show the redshift distribution of the galaxy catalog.

Fig. 3. Histograms of galaxy-catalog redshifts (right) and of corresponding luminosity distances in the three MG scenarios considered in this work (left). MG0.6 is the sample with Ξ0 = 0.6, GR with Ξ0 = 1, and MG1.8 with Ξ0 = 1.8.

The sample of GW events was created by populating the galaxy catalog with CBC events. The sky position and redshift of each event are determined by the host galaxy, while the luminosity distance was obtained by assuming fiducial cosmological, and MG propagation models. We chose the same cosmological model adopted in the MICE simulation: a flat ΛCDM universe with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m, 0 = 0.25. We parameterized MG propagation with two parameters (Ξ₀, n) as in Eq. (2). Binary masses in the source frame were drawn from a mass distribution factorized as

$$\begin{array}{c}\hfill p(m_1,m_2\mid {\mathit{\lambda }}_m)=p(m_1\mid {\mathit{\lambda }}_m)p(m_2\mid m_1,{\mathit{\lambda }}_m).\end{array}$$(20)

The primary mass follows a power law with spectral index −α, truncated in the [m_low, m_high] range, with an additional Gaussian peak centered at μ_g with a width of σ_g, weighted by λ_g. The lower edge of the power law is smoothed using a parameter δ_m. The secondary mass follows a smoothed power law with spectral index β, truncated in the [m_low, m₁] range. The chosen mass model, fully described in Appendix B of Abbott et al. (2021c), is consistent with the latest GWTC-3 constraints (Abbott et al. 2023a). For the merger-rate evolution, we used the Madau-Dickinson parameterization (Madau & Dickinson 2014; Fishbach & Holz 2017):

$$\begin{array}{c}\hfill \psi (z;{\mathit{\lambda }}_r)\propto \frac{{(1+z)}^{\gamma }}{1+{\left(\frac{1+z}{1+z_p}\right)}^{\gamma +\kappa }}·\end{array}$$(21)

For the cosmological, rate, and mass hyperparameters a single fiducial value was chosen. For the MG hyperparameters, three combinations are explored, resulting in three different GW populations. One case, Ξ₀ = 1, corresponds to GR with no modified GW propagation. The other cases, Ξ₀ = 0.6 and Ξ₀ = 1.8, are at the boundaries of the 1-σ constraints obtained by Mancarella et al. (2022) using O3 data and considering only features of the population models. The value of n is set to 2.7 in all scenarios. Starting from the same redshift distribution, the corresponding luminosity distance distributions of the three populations are expected to span different ranges due to the different values of Ξ₀, affecting both the number and the distributions of detected events. These differences can be inspected in the right panel of Fig. 3, where we show the luminosity distance distributions of the galaxy catalog in the three different MG scenarios.

Table 1 summarizes the hyperparameters describing the GW populations, their fiducial values, and the prior used in the following MCMC analyses.

3.2. GW detection and PE generation

The sample of GW events was analyzed using GWFAST (Iacovelli et al. 2022a,b) to identify detectable events. We assume the CBC event to be quasi-circular, non-precessing BBH system. Each CBC waveform is characterized by 15 detector-frame parameters:

$$\begin{array}{c}\hfill {\mathit{\theta }}^{\mathrm{d}}=\{{\mathcal{M}}_c,\eta ,d_L,\theta ,\varphi ,\iota ,{\chi }_1^z,{\chi }_2^z,\varphi ,t_c,{\mathrm{\Phi }}_c\},\end{array}$$(22)

where ℳ_c is the redshifted chirp mass, η is the symmetric mass ratio, d_L is the binary luminosity distance, (θ, ϕ) are the sky position angles, ι is the inclination angle between the orbital angular momentum and the line of sight, χ_1, 2^z are the spin projections along the orbital angular momentum, ψ is the polarization angle, t_c is the coalescence time, and Φ_c is the phase at coalescence. While the first five parameters were generated based on the galaxy catalog and the population models as explained in the previous subsection, the remaining parameters were drawn from specific distributions: the spin components are uniformly distributed in [ − 1, 1], the inclination angle is uniform in cosι over [0, π], the polarization angle and coalescence phase are uniformly distributed in [0, π] and [0, 2π], respectively, and the coalescence time – expressed in units of fraction of a day – is uniform in [0, 1].

The waveform signal, simulated using the IMRPhenomHM (London et al. 2018) approximant, is injected into a noise realization of each detector. We considered a network configuration that includes the two LIGO interferometers in the USA (Aasi et al. 2015), the Virgo interferometer in Italy (Acernese et al. 2015), the KAGRA interferometer in Japan (Aso et al. 2013), and the planned LIGO interferometer in India (LIGO 2011). We assumed the publicly available² sensitivity curves representative of the O5 observing run (Abbott et al. 2016b), with 100% duty cycle. We then analyzed the injected signals with GWFAST, estimating the network’s match-filtered signal-to-noise ratio (S/N) and computing the Fisher-information-matrix (FIM) and its inverse. The individual GW event likelihood for the detector-frame parameters is approximated as a multivariate Gaussian, with the covariance matrix given by the inverse FIM. In Borghi et al. (2024), we checked that the Fisher matrix approach was a good approximation for high S/N events, such as the ones considered in the following analyses. The PE samples for θ^d are drawn from this approximated likelihood using the emcee sampler. The priors used are uniform in the [0, 10⁵] M_⊙, [0, 1/4] range, and in [0, 10⁵] Gpc for ℳ_c, η, and d_L, respectively, while all other waveform parameters were bounded in the same physical ranges from which they were drawn. The prior also includes the Jacobian of the transformation (ℳ_c, η)→(m₁^d, m₂^d) to ensure that the binary masses PE samples are uniformly distributed.

To create the injection set, we adopted a similar process, excluding the unnecessary PE generation step. We used the same N_inj = 2 × 10⁷ injections as in the O5-like mock catalog described in Borghi et al. (2024).

3.3. GW catalogs selection and properties

Based on the population models described in Eq. (3.1) and assuming a local BBH merger rate density of R₀ = 17 Gpc⁻³ yr⁻¹ (Abbott et al. 2023a), the expected number of events per year out to z = 1.3 (the upper limit of our galaxy catalog) is:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{\mathrm{cbc}}}{\mathrm{d}{t}_{\mathrm{d}}}=R_0{\displaystyle {\int }_0^{1.3}}\mathrm{d}z\frac{\psi (z;{\mathit{\lambda }}_r)}{1+z}\frac{\mathrm{d}{V}_c}{\mathrm{d}z}(z;{\mathit{\lambda }}_c)\approx 1.4\times {10}^4{\mathrm{yr}}^{-1}.\end{array}$$(23)

We note that this rate is strongly affected by the chosen value of R₀, which is still largely unbounded from O3 data $R_0={17}_{-6.7}^{10}{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}$. From this total, the number of detected events with S/N > 8 is approximately 6200, 3100, and 1100 in the MG0.6, GR, and MG1.8 scenarios, respectively. This difference can be understood from the luminosity distance distributions of the three GW samples in Fig. 3, as higher (lower) Ξ₀ values map the true redshifts into higher (lower) luminosity distances, affecting the detection rates. Similarly, the number of events with S/N > 20 per year is approximately 915, 340, and 100 in the three scenarios.

We selected 300 events at S/N > 20 from each GW catalog. This corresponds to roughly one year of observation assuming that GR holds, or three years (4.5 months) for Ξ₀ = 1.8 (Ξ₀ = 0.6). Keeping the same number of events enables a fair comparison of the derived constraints, which strongly depend on the number of events. On the other hand, since different values of Ξ₀ have an impact on the number of GW events detected per year, we also analyzed results within a fixed time frame using two additional catalogs: 900 sources for MG0.6 and 100 for MG1.8. In Section 4, we discuss the impact on the results of this different scenario.

The properties of the three catalogs are illustrated in Fig. 4. While the primary mass distribution appears similar across the three cases, the redshift histogram of the MG1.8 catalog shows more events at lower redshifts. Again, this effect can be understood since high-redshift GW events are less likely to be detected as they correspond to very high-luminosity distances in this scenario. For the opposite reason, the MG0.6 catalog extends to higher redshifts. Since the cut in S/N is the same for all three catalogs, the luminosity distance and localization area measurement uncertainties are similar. The events in the MG1.8 catalog typically have fewer galaxies in the localization volume, which is expected due to the lower number of galaxies at lower redshifts. In particular, this catalog contains seven “golden sirens”, defined here as BBHs with 100 or fewer galaxies within their localization volume. In comparison, the MG0.6 and GR catalogs have one and three golden sirens, respectively, since their redshift distributions are shifted to higher values.

Fig. 4. Properties of three GW mock catalogs. From left to right: primary mass distributions, redshift distributions, luminosity distance, localization area uncertainties, and histograms of the number of galaxies in the localization volume.

4. Results

In this section, we discuss the results obtained. For each GW catalog, we considered both spectroscopic (spec-z) and photometric uncertainties on galaxy redshifts:

$$\begin{array}{c}\hfill {\stackrel{\sim }{\sigma }}_{z,g}=\{\begin{array}{cc}0.001(1+z),\hfill & \mathtt{spec}-\mathtt{z};\hfill \\ 0.05(1+z),\hfill & \mathtt{photo}-\mathtt{z}.\hfill \end{array}\end{array}$$(24)

Spectroscopic galaxy catalogs can be obtained by expanding the currently available catalog (GLADE+ Dálya et al. 2022) with future data. For example, the all-sky ESA Euclid survey (Laureijs et al. 2011) will measure spectroscopic redshift in the 0.9 < z < 1.8 range with an accuracy level of ${\stackrel{\sim }{\sigma }}_{z,g}/(1+z)\lesssim 0.001$; DESI (DESI Collaboration 2016) was planned to observe about one-third of the sky and cover the redshift range of 0.4 < z < 2.1; the WST will map roughly half of the sky up to a target redshift of 1.5 (Mainieri et al. 2024). Photometric redshifts are also available with ongoing surveys. For instance, the DES survey reached ${\stackrel{\sim }{\sigma }}_{z,g}\sim 0.01$ (Myles et al. 2021) over a smaller area, with the potential to improve to ${\stackrel{\sim }{\sigma }}_{z,g}\sim 0.007$ using advanced techniques (Buchs et al. 2019). Euclid and the upcoming Rubin observatory (Ivezić et al. 2019) are expected to provide photometric redshifts with an accuracy of ${\stackrel{\sim }{\sigma }}_{z,g}\sim 0.05(1+z)$ (Desprez et al. 2020; Schirmer et al. 2022).

For each GW and galaxy catalog (with different redshift assumptions), we sampled the posterior with an MCMC approach, exploring the following configurations:

1. inference on all hyperparameters listed in Table 1;

2. inference only on MG and population hyperparameters, fixing H₀ to its fiducial value;

3. inference on H₀ and population hyperparameters, fixing MG ones to their fiducial values.

In all cases, Ω_m, 0 is fixed. For each GW catalog, this results in six MCMC fits for a total of 18 runs. The priors used are summarized in Table 1. We used the emcee sampler with 100 walkers and evolved the chains until the number of samples was at least 50 times larger than the integrated autocorrelation time for all the hyperparameters. The “many-1D” kernel is employed, enabling each fit to converge in about 18 hours on a single GPU. This extensive series of 18 tests, which cumulatively represent a population of about 5000 events, would not have been possible with CHIMERA 1.0, which was already at its limit with only 100 sourced events. Using CHIMERA 1.0 for this work would have taken approximately 270 CPU days or 18 GPU months (see Fig. 2).

4.1. Hyperparameter constraints

The left panel of Fig. 5 shows the 1D marginalized posterior distributions obtained in the first MCMC configuration using spectroscopic galaxy redshifts for the three GW catalogs. The constraints are shown for cosmological, MG, and selected population hyperparameters. For the latter, we focused on the position and width of the Gaussian mass peak and the low-z slope of the rate. Indeed, these are the hyperparameters that mostly correlate with the cosmological and MG ones. The right panel of Fig. 5 shows the corresponding results using photometric galaxy redshifts. The colored areas represent the 68% confidence levels (C.L.s) around the median of the distributions. Table 2 summarizes the median, the 68% C.L., and the percentage precision for cosmological and MG hyperparameters across all MCMC configurations. The percentage precisions on H₀ and Ξ₀ across the three catalogs and the different MCMC configurations are also shown in Fig. 6. To visualize how constraints on MG and cosmological hyperparameters vary across the three GW catalogs and between photometric and spectroscopic redshifts, we project them on the luminosity distance-redshift relation (2), as shown in Fig. 7.

Fig. 5. Marginalized posterior distributions for selected hyperparameters. These constraints are obtained in the first MCMC configuration (inference on all hyperparameters) and assuming a spectroscopic (left panel) or photometric (right panel) galaxy catalog. The colored areas under each posterior represent the 68% C.L. The dotted lines indicate the hyperparameter fiducial values.

Fig. 6. Percentage precision on H0 and Ξ0 obtained with three GW catalogs across the various MCMC configurations. Dark (light) gray band represents the 1% (10%) error range.

Fig. 7. Projected constraints on dLgw − z relation (2) obtained in first MCMC configuration. The darker (lighter) contours represent the 68% (95%) C.L.

We now discuss the results with spectroscopic galaxy redshifts. H₀ and Ξ₀ were recovered without bias at 68% C.L. across all MCMC configurations. In the second MCMC configuration, Ξ₀ was constrained with a precision of 22%, 7.5%, and 10% for Ξ₀ = 0.6, 1, and 1.8, respectively. This allows the three scenarios to be distinguished at 68% C.L., as shown in the right panel of Fig. 5. Notably, fixing H₀ does not improve constraints on the MG hyperparameters (see Table 2).

Fig. 8. 2D contours between cosmological, MG, and selected population hyperparameters across the different GW and galaxy catalogs. The 68% and 95% confidence level regions are indicated by dark and light colors, respectively. The dotted lines represent the hyperparameter fiducial values.

The precision on H₀ depends on whether MG hyperparameters are marginalized or fixed. When the MG hyperparameters were fixed to their fiducial values, H₀ was recovered with a precision on the order of or slightly better than 1%. The best constraint, 0.72%, was obtained in the MG1.8 case, likely due to its higher number of golden sirens compared to the other catalogs. In fact, when the MG hyperparameters are not inferred, the precision on H₀ improves with the number of golden sirens in the catalog. However, when MG hyperparameters are marginalized, the precision on H₀ degrades by factors of ∼2.1, ∼3.9, and ∼7.2 for the MG0.6, GR, and MG1.8 cases, respectively. In this scenario, the best result is found with the MG0.6 catalog, suggesting that the H₀ precision is no longer strongly influenced by the number of golden sirens.

The posterior distributions for n remain nearly flat, even when H₀ is not marginalized.

With regard to the results with photometric galaxy redshifts, switching from spectroscopic to photometric redshifts does not introduce biases, but significantly weakens the constraints on H₀ and Ξ₀. Specifically, the precision on Ξ₀ degrades by a factor of ∼1.4, ∼5.4, and ∼3.5 for the MG0.6, GR, and MG1.8 cases, respectively. Fixing H₀ mitigates this degradation, and Ξ₀ is recovered with a precision that is on average only 1.2 times worse than the spectroscopic case with fixed H₀. The precision on H₀ degraded by factors of ∼5 − 7 when the MG hyperparameters were fixed, and by ∼4 − 10 when Ξ₀ and n were marginalized. As in the spectroscopic case, n remains unconstrained in the photometric scenario.

Concerning constraints for a fixed time frame, we tested how constraints change when considering a fixed time frame of a one-year of observation. To do this, we used the additional catalogs mentioned above of 900 and 100 sources for the MG1.8 and MG0.6 cases. We performed a test considering spectroscopic galaxies and inferring both MG and cosmological hyperparameters. We find that in the MG0.6 case, the constraint on Ξ₀ improves from 22% to 17%, while the precision on H₀ remains approximately the same. In the MG1.8 case, instead, the lower number of events per year degrades the precision of Ξ₀ from 10% to 14% and of H₀ from 5.2% to 8.5%.

With regard to population hyperparameters, the Gaussian peak position and width of the primary mass distribution are recovered within the 68% CL across all combinations of GW catalogs and redshift error assumptions. The redshift error choice has only a marginal effect on the posteriors, and fixing either MG or cosmological hyperparameters has no impact on these hyperparameters. This implies that population hyperparameter constraints are mainly driven by the GW data rather than the galaxy catalog used.

Similar results were obtained for the merger rate parameter γ, though the fiducial value lies slightly outside the 68% C.L. for MG0.6 and GR. Several factors may contribute to this: the MG1.8 catalog contains a higher density of events at low redshifts, potentially enhancing its constraining power; alternatively, the problem could be due to an insufficient parameter-space sampling in the current injection set. Nevertheless, since both H₀ and Ξ₀ are recovered without bias and exhibit a much weaker correlation with γ than with each other (see next section and Fig. 9), this deviation does not affect our conclusions. A more detailed investigation of this issue is left for future work.

Fig. 9. Correlation among selected hyperparameters in first MCMC configuration. The strength and direction of the correlation were quantified with a Spearman coefficient. In each panel, the upper corner plot shows the results in the photo-z case, while the lower triangle shows the results in the spec-z case. The left panel shows results from the MG0.6 GW catalog, the center panel from the GR catalog, and the right panel from the MG1.8 GW catalog.

4.2. Hyperparameter correlations

In Fig. 8, we show the correlations between cosmological, MG, and selected population hyperparameters across different catalogs and redshift error assumptions. For the population hyperparameters, we focused on the mass peak position, μ_g, and on the first slope of the merger rate, γ, as these are the most correlated with H₀ and Ξ₀. To further quantify these correlations, we calculated the Spearman rank correlation coefficient (ρ) (Zwillinger & Kokoska 1999). This metric, which ranges from −1 to 1, quantifies the strength and direction of the correlation: a positive value indicates correlated variables, a negative value indicates a monotonically decreasing relation, and zero represents two non-correlated variables. Fig. 9 shows Spearman coefficients obtained when inferring all hyperparameters for the three GW catalogs. The lower (upper) triangle of each panel displays the results obtained using spectroscopic (photometric) redshifts.

In the photometric case, all catalogs exhibit a strong correlation between Ξ₀ and H₀. When spectroscopic redshifts are used, the degeneracy between these two hyperparameters is significantly reduced – except in the GR case. Furthermore, the constraints become much tighter due to the additional information provided by the spectroscopic galaxy catalog.

In the MG0.6 spec-z, n is strongly correlated with Ξ₀ and anticorrelated with H₀. These correlations reverse sign in the MG1.8 spec-z case. In contrast, for the GR spec-z scenario, n shows no significant correlation with either H₀ or Ξ₀. We also note that in all photo-z cases, the Spearman coefficients of n with H₀ and Ξ₀ are smaller than the in the respective spec-z cases; however this is due to the much weakerconstraints on n.

When using photometric redshifts, we find that Ξ₀ is slightly correlated with the mass peak μ_g – particularly in the GR case – and strongly correlated with γ. However, the Ξ₀-γ correlation is notably weaker than when the galaxy catalog is not included (see, e.g., the corner plots 5 and 6 in Mancarella et al. 2022). Remarkably, we find that using spectroscopic redshifts breaks these correlations.

5. Conclusions

In this work, we presented an enhanced version of CHIMERA, a code designed for Bayesian inference of cosmological, modified-gravity, and population hyperparameters using standard sirens and galaxy catalogs. This second release addresses the computational bottlenecks of the first version by including three distinct KDE algorithms used to compute the likelihood. We introduced these algorithms, detailing their differences in terms of results and computational efficiency, and validated their results against one another. As a next step, we plan to extend this validation test to real data and/or larger datasets to assess possible systematic effects. This will be one of the main goals of a future blinded-mock-data challenge between CHIMERA 2.0 and similar pipelines.

Using CHIMERA 2.0, we forecasted constraints for an O5-like detector network in an optimistic scenario. Three BBH populations have been studied: one assuming no modified GW propagation, one where the effective luminosity distance is greater than the electromagnetic one (Ξ₀ = 1.8) and one where it is smaller (Ξ₀ = 0.6). Population hyperparameters are jointly inferred with MG and cosmological ones. We studied both cases in which the galaxies have spectroscopic and photometric redshifts. The fiducial MG and cosmological hyperparameters are recovered at 68% C.L. in all scenarios. With spectroscopic redshift, the fiducial values of Ξ₀ can be distinguished across all scenarios, regardless of whether H₀ is marginalized or fixed. However, the H₀ precision, which is always below 1% when the MG hyperparameters are fixed, degrades by factors of ∼2.1, ∼3.9, and ∼7.2 when the MG hyperparameters are marginalized. Notably, the marginalization over MG hyperparameters strongly weakens the dependence of the H₀ precision on the number of golden sirens. In the photometric case, the three Ξ₀ values can only be distinguished when H₀ is not inferred, and become indistinguishable when H₀ is marginalized. In the photometric case, the constraints on H₀ degrade on average by a factor of 3.4 (1.2) when the MG hyperparameters are marginalized (fixed) compared to the corresponding spectroscopic case. In addition to this, using spectroscopic redshifts significantly reduces the correlation among cosmological, MG, and population hyperparameters. This again highlights the importance of future spectroscopic surveys, such as WST (Mainieri et al. 2024), in terms of fully exploiting GWs as standard sirens for probing cosmology and modified GW propagation. Lastly, the posterior distributions for n remain nearly flat, and fixing H₀ does not improve the constraints on this parameter, both in the spectroscopic and photometric case.

We emphasize that these results are based on several simplifying assumptions. First, we assume a complete galaxy catalog that includes only the most massive galaxies. In real galaxy surveys, however, galaxy-selection functions are more complex and must be properly accounted for in the cosmological analysis. At the same time, the completeness correction must be properly modeled, reflecting the probability of a galaxy hosting a merger event based on its astrophysical properties, rather than relying on uninformative assumptions such as uniform comoving volume completion. A more realistic and physically motivated framework that addresses these complexities will be developed in a future study. While our current implementation uses this simplified framework, the results presented in this paper highlight the need for spectroscopic galaxy surveys and multiple GW detectors able to precisely localize GW sources in standard siren cosmology and GR testing. Lastly, we plan to expand these forecasts to next-generation interferometers and prove how constraints will improve. Although CHIMERA 2.0 can handle approximately 50 times more events than CHIMERA 1.0, handling about 10⁵ events would require additional computational power. This could be achieved by parallelizing the likelihood calculation across multiple GPUs using JAX functionalities or the Message Passage Interface, which is a topic for future work.

Acknowledgments

We thank Michele Mancarella for useful discussions and comments. We acknowledge the ICSC for awarding this project access to the EuroHPC supercomputer LEONARDO, hosted by CINECA (Italy). This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation. MT acknowledges the funding from the European Union – NextGenerationEU, in the framework of the HPC project – “National Center for HPC, Big Data and Quantum Computing” (PNRR – M4C2 – I1.4 – CN00000013 – CUP J33C22001170001). MM acknowledges the financial contribution from the grant PRIN-MUR 2022 2022NY2ZRS 001 “Optimizing the extraction of cosmological information from Large Scale Structure analysis in view of the next large spectroscopic surveys” supported by NextGenerationEU. MM and NB acknowledge the financial contribution from the grant ASI n. 2024-10-HH.0 “Attività scientifiche per la missione Euclid – fase E”. We acknowledge the use of the following software: NumPy (Harris et al. 2020), JAX (Bradbury et al. 2018), equinox (Kidger & Garcia 2021), matplotlib (Hunter 2007), seaborn (Waskom 2021), arviz (Kumar et al. 2019), emcee (Foreman-Mackey et al. 2013), ChainConsumer (Hinton 2016), GWFAST (Iacovelli et al. 2022b), CHIMERA 1.0 (Borghi et al. 2024).

References

Appendix A: Kernels comparison and validation

In Fig. A.1 we compare the MCMC results obtained using three different kernels: the one that implements CHIMERA 1.0, the "many-1D" (18), and the "single-1D" (19). The results are obtained using the “O5-like” catalog of Borghi et al. (2024) and spectroscopic galaxy redshifts. Overall, the posteriors obtained in the various cases are in excellent agreement.

Fig. A.1. Comparison between the MCMC results obtained using the kernels implemented in CHIMERA 2.0 and the one implemented in CHIMERA 1.0. The contours represent the 68% and 95% C.L. The dotted lines indicate the hyperparameter fiducial values.

All Tables

All Figures

	Fig. 1. Top: Pixelization of 90% localization area of a mock GW event, with PE samples (marked with a cross) colored according to the pixel they fall into. Bottom: Different GW kernels implemented in CHIMERA 2.0. Each line represents the GW kernel in a pixel of the above mock GW event. Kernels are computed at single hyperparameter space points.
In the text

Fig. 2. Top: Likelihood evaluation times for catalogs with different numbers of events (each with 5000 PE samples) at a single hyperparameter space point, compared across different kernel algorithms. Solid curves show times on four Intel CPUs (2.0 GHz clock frequency), while dashed curves represent times on a single Nvidia A100 GPU. Bottom: Run times for a full affine-invariant MCMC fit, using the emcee sampler with 50 walkers, on 25 CPUs (left panel) or one GPU (right panel).

In the text

	Fig. 3. Histograms of galaxy-catalog redshifts (right) and of corresponding luminosity distances in the three MG scenarios considered in this work (left). MG0.6 is the sample with Ξ0 = 0.6, GR with Ξ0 = 1, and MG1.8 with Ξ0 = 1.8.
In the text

	Fig. 4. Properties of three GW mock catalogs. From left to right: primary mass distributions, redshift distributions, luminosity distance, localization area uncertainties, and histograms of the number of galaxies in the localization volume.
In the text

	Fig. 5. Marginalized posterior distributions for selected hyperparameters. These constraints are obtained in the first MCMC configuration (inference on all hyperparameters) and assuming a spectroscopic (left panel) or photometric (right panel) galaxy catalog. The colored areas under each posterior represent the 68% C.L. The dotted lines indicate the hyperparameter fiducial values.
In the text

	Fig. 6. Percentage precision on H0 and Ξ0 obtained with three GW catalogs across the various MCMC configurations. Dark (light) gray band represents the 1% (10%) error range.
In the text

	Fig. 7. Projected constraints on dLgw − z relation (2) obtained in first MCMC configuration. The darker (lighter) contours represent the 68% (95%) C.L.
In the text

	Fig. 8. 2D contours between cosmological, MG, and selected population hyperparameters across the different GW and galaxy catalogs. The 68% and 95% confidence level regions are indicated by dark and light colors, respectively. The dotted lines represent the hyperparameter fiducial values.
In the text

Fig. 9. Correlation among selected hyperparameters in first MCMC configuration. The strength and direction of the correlation were quantified with a Spearman coefficient. In each panel, the upper corner plot shows the results in the photo-z case, while the lower triangle shows the results in the spec-z case. The left panel shows results from the MG0.6 GW catalog, the center panel from the GR catalog, and the right panel from the MG1.8 GW catalog.

In the text

	Fig. A.1. Comparison between the MCMC results obtained using the kernels implemented in CHIMERA 2.0 and the one implemented in CHIMERA 1.0. The contours represent the 68% and 95% C.L. The dotted lines indicate the hyperparameter fiducial values.
In the text