Non-Parametric Reconstruction of the Hubble Parameter from the Fourth Gravitational Wave Transient Catalog and DESI Baryonic Acoustic Oscillations

IOP Publishing Journal vv (yyyy) aaaaaa Author et al Classical Quantum Gravity ARTICLE Crossmark RECEIVED dd Month yyyy REVISED dd Month yyyy Non-Parametric Reconstruction of the Hubble Parameter from the Fourth Gravitational Wave Transient Catalog and DESI Baryonic Acoustic Oscillations Grégoire Pierra1,∗, Alberto Colombo1,2and Simone Mastrogiovanni1 1INFN, Sezione di Roma, 1-00185 Roma, Italy 2INAF – Osservatorio Astronomico di Brera, via Emilio Bianchi 46, I-23807 Merate (LC), Italy ∗Author to whom any correspondence should be addressed. E-mail: [email protected] Keywords: cosmic expansion history, gravitational-wave, BAO, cosmology, spectral siren, non-parametric inference Abstract The release of the fourth Gravitational-Wave Transient Catalog (GWTC-4.0) by the LIGO–Virgo–KAGRA collaboration includes more than 200 compact binary coalescence (CBC) candidates that can be used to probe the cosmic expansion. The population of merging binary black holes has been used so far to provide a constraint on the Hubble constant and dark matter fraction under the hypothesis of a flat-Λ-Cold-Dark-Matter Universe. In this work, we provide the first non-parametric constrain on the Hubble parameter from 137 dark sirens reported in GWTC-4.0. We employ the relation between detector and source frame masses for detected GW signals, to obtain a statistical redshift evaluation for the population of binary black holes (BBHs). We model the Hubble parameter as a non-parametric autoregressive process in terms of the scale factor, using splines. In addition, we introduce two novel features: the use of anchor points for H(z)derived from an external probe — here, Baryon Acoustic Oscillations (BAOs) — and a constraining power coefficient that quantifies where the inference is most data-driven by GW detections. We highlight three key findings: (i) using GWs alone, the Hubble parameter determination is the most GW-data-driven around redshift z= 0.44, yielding to H(0.44) = 92.3+29.9 −36.6kms−1Mpc−1. Its value at z= 0, the Hubble constant, is therefore less constrained by the GW data. (ii) The Hubble parameter inferred from analyses assuming a flat-ΛCDM cosmological model is strongly affected by the cosmological model assumption. (iii) Introducing an anchor point for H(z)enhances the inferred constraints and provides a clear visualization of the redshift range where GWs contribute most to the constraining power. 1 Introduction Measuring the cosmic expansion history of the Universe is a central goal of modern cosmology, with a particular focus on determining its local expansion rate (H0), the Hubble constant. This question has attracted considerable attention due to the growing discrepancy between early–universe estimates of H0inferred from the cosmic microwave background (CMB) [1] and late–universe measurements based on standardized candles such as type Ia supernovae [2]. This discrepancy, widely referred to as the “Hubble tension”, has now surpassed the 5σlevel. Gravitational waves (GWs) from compact binary coalescences (CBCs) have recently emerged as an independent and powerful probe of the Hubble constant, complementing both CMB and standard-candle-based methods. The concept of using GWs as standard sirens, by analogy with standard candles in cosmology, was first proposed by Schutz in 1986 [3] and has since been extended in numerous studies when an electromagnetic counterpart is present [4–7], and when not [3,8–21]. GW observations provide a direct measurement of the luminosity distance to the source, independent of the traditional distance ladder, and when paired with a redshift measurement, enable constraints on cosmological parameters such as the Hubble constant. In this paper, we focus 1 arXiv:2511.11795v1 [astro-ph.CO] 14 Nov 2025 IOP Publishing Journal vv (yyyy) aaaaaa Author et al on the “spectral siren” method [13,22], a technique able to exploit the difference between detector and source mass and obtain a redshift estimation for the population of GW sources. While determining H0is important, an even more compelling objective is to chart the expansion history of the Universe through the Hubble parameter, H(z). Measuring H(z)as a function of redshift allows us to constrain cosmological models, investigate the evolution of the expansion, and potentially probe the nature of dark energy and dark matter. GWs from binary black hole observations are excellent sources for this task as they are observed up to redshift 2 with current detectors sensitivities. Current cosmological expansion measurements from GW sources typically adopt flat-ΛCDM framework [21,23–28], and thus the reconstructed behavior of H(z)is implied by the model assumption and only constraints on H0or Ωmcan be reported. With the recent release of the largest GW catalog to date by the LIGO–Virgo–KAGRA (LVK) collaboration [29–31], the Gravitational-Wave Transient Catalog 4.0 (GWTC-4.0) [31–33], which contains a total of 218 candidate detections from CBCs, the potential of standard sirens for cosmological studies can now be explored in unprecedented detail. We present the first non-parametric spectral siren inference of the Hubble constant, and more generally of the Hubble parameter H(z), using GW observations. Specifically, we analyze 137 BBHs with a false-alarm rate (FAR) below 0.25 yr−1, as reported in GWTC-4.0 [34], modeling the Hubble parameter as an autoregressive process with splines. We further assess the extent to which GW data themselves constrain the expansion history across redshift, compared with the constraints induced by model assumptions in parametric approaches. Our proposed model further allows for an agnostic combination of GW sources detected across a broad luminosity-distance range with independent constraints on the Hubble parameter H(z)from other cosmological probes. In this work, we specifically investigate the effect of incorporating external H(z)measurements, for instance those provided by Baryonic Acoustic Oscillations (BAOs) from DESI [35]. Some studies have also investigated the use of BAO measurements in the context of testing modified gravity theories, although used in a different framework with a different application from the one explored in this work [36–38]. The paper is organized as follows. In Section 2, we present our new approach for constructing a non-parametric model of the Hubble parameter using splines and discuss its connection to the equation of state of the Universe. Section 3introduces the Bayesian framework used to perform the spectral siren inference, along with the population models and pipeline employed for the analyses. Section 4summarizes our main findings, including the non-parametric estimation of the Hubble parameter from GWs, with and without anchors from BAOs, and comparison with the widely used flat-ΛCDM cosmological model. It also quantifies the effect of using parameterized cosmological models for spectral sirens and identifies when the inference of H(z)is genuinely driven by the GW data. Finally, we conclude in Section 5. 2 A non-parametric cosmological model The spectral siren method estimates the cosmological redshifts of GW sources by exploiting the intrinsic relation between the detector mass mdand the source mass ms: md= (1 + z)ms,(1) where zis the cosmological redshift of the source, and mdis commonly referred to as the redshifted mass. In practice, however, the redshift of a GW source is not directly measured from the signal, unlike the luminosity distance, but is instead inferred statistically. This is achieved by fitting the distribution of BBHs in the luminosity distance–detector mass space, assuming a parametric source mass model and using the cosmological evolution encoded in the Friedmann equations. Figure 1illustrates the procedure underlying the spectral siren inference method. It shows the evolution of two mass spectra in the luminosity distance–detector-frame mass plane, assuming a flat-ΛCDM cosmological model with two different choices for the Hubble constant and matter-density parameters. For nearby sources (dL≈0 Mpc), i.e., GW signals emitted close to the observer, the detector mass spectra coincide for both cosmologies. This is expected because, according to Eq. 1, the source mass equals the detector mass at vanishing redshift. As sources are pushed to larger distances, the detector mass spectrum evolves, reflecting its dependence on the cosmological redshift z. This evolution occurs because the source mass distribution is assumed to be redshift-independent, while the mapping between detector and source masses depends on the cosmological parameters. For illustration, the shift of the Gaussian peak initially located at 35 M⊙ for dL= 0 Mpc highlights how the evolution of the mass spectrum varies with the adopted cosmological model and parameter values. 2 IOP Publishing Journal vv (yyyy) aaaaaa Author et al 0 20 40 60 80 100 120 140 Detector mass md[M] 0 1.25 2.50 3.75 5.00 Luminosity distance dL[Gpc] H0= 70 km s Mpc ; Ωm,0= 0.3 H0= 200 km s Mpc ; Ωm,0= 0.8 Figure 1. Graphical representation of the evolution of the mass spectrum in the detector frame as observed at increasing luminosity distances, assuming different values of the Hubble constant (H0) and the matter-density parameter (Ωm,0) within a flat-ΛCDM cosmological model. The figure illustrates how changes in cosmological parameters affect the observed shape of mass features at larger distances. The evolution of the mass spectrum for a given cosmology is governed by the Friedmann equations, and in particular by the redshift dependence of the Hubble parameter H(z), since md= (1 + z)ms= [1 + f(dL;H(z))] ms,(2) where f(·)denotes the mapping between luminosity distance and cosmological redshift for a specified cosmological model, here assumed to be a flat-ΛCDM Universe. Throughout the remainder of this study, we adopt this flat-ΛCDM framework as the fiducial model, or stated otherwise when the non-parametric approach is employed. In the low-redshift regime, the function f(·)can be approximated as f(dL;H(z=0)) ≈dLH0 c,(3) where cis the speed of light and H0denotes the present-day value of the Hubble parameter, the Hubble constant. A direct consequence within flat-ΛCDM cosmological models is that higher values of H0lead to a more pronounced evolution of the detector mass spectrum, as illustrated in Fig. 1. In the low-redshift limit, however, the matter-density parameter Ωm,0plays no significant role, as indicated by Eq. 3. In spectral siren cosmology, the constraining power of this approach arises from the fact that, for a given luminosity distance, different cosmological models predict different redshifts through Eq. 2, resulting in distinct detector-frame mass distributions and evolutionary trends. In this study, we aim to develop a non-parametric method to infer the cosmic expansion history solely from GW observations. Accordingly, it is natural to represent the cosmological model in a non-parametric form by mapping selected redshift values to corresponding luminosity distance nodes or bins. Since GW signals are observed in the luminosity distance–detector mass parameter space, this framework allows us to identify which GW candidates contribute to the estimation of the Hubble parameter, now expressed as a function of luminosity distance, H(dL). These non-parametric measurements can then be compared a posteriori with the predictions of standard cosmological models, such as the flat-ΛCDM Universe. However, adopting a non-parametric model of luminosity distance as a function of redshift introduces significant challenges. First, the spectral siren method requires the Hubble parameter to estimate the rate density of CBC mergers in the Universe. A non-parametric representation of z(dL)can, in some cases, yield negative values of the Hubble parameter, even if z(dL)is 3 IOP Publishing Journal vv (yyyy) aaaaaa Author et al constrained to be monotonic. Such behavior would correspond to a Universe that contracts rather than expands, which is physically inconsistent with current observations. Moreover, certain non-parametric functions, such as splines [39–42], one of the simplest and most commonly used non-parametric representations, can produce discontinuities in H(z)at the redshift nodes. Second, depending on how the function z(dL)is constructed, the method may explore cosmologies with extreme or otherwise unphysical values of H(z), for instance in cases where z(dL)evolves very slowly with distance. Given the challenges discussed above, we construct our non-parametric model directly for the Hubble parameter, H(z), as a function of redshift, using a spline parameterization in natural logarithmic space. The cosmological redshift range is divided into Nnodes redshift nodes, zi, and the natural logarithm of the Hubble parameter is then defined such that ln H(z) = Nnodes−1 X i=1 ailn 1+z 1+zi+biΘ(z;zi, zi+1),(4) where Θ(·)is a box function ensuring no contribution to the Hubble parameter from intervals outside [zi, zi+1]. From Eq. 4, the coefficients aiand bican be expressed as ai= ln H(zi+1) H(zi)/ln 1+zi+1 1+zi(5) bi= ln H(zi).(6) This spline-based parameterization offers several advantages. First, the Hubble parameter constructed in this way is a continuous function, in contrast to the discontinuities encountered in other non-parametric approaches. Second, and more importantly from a physical perspective, the Hubble parameter now scales locally as H(z)∝(1+z)α=aα, where ais the scale factor thus providing a direct connection to the equation of state (EOS) of the Universe and its energetic content. Indeed, under the hypothesis that a single energy component dominates within a redshift bin, the Hubble parameter can be rewritten as H(z)=H(zi)1+z 1+zi3(1+wi) =H(ai)hai ai3(1+wi),(7) where, for example, wi= 0 corresponds to a Universe dominated by cold dark matter in the redshift bin zi. Using the parametrization introduced in Eq. 4, the comoving volume element between redshifts ziand zi+1, for zi< z < zi+1, can be written as ∆dcom,i(z) =    c(1+zi) H(zi)(1−ai)(1+z)1−aih1+zi+1 1+z H(z) H(zi+1)−1iif ai= 1, 1+zi H(zi)ln h1+zi+1 1+ziif ai= 1 .(8) This expression provides a convenient way to approximate the comoving distance by discretizing the redshift range into intervals defined by the nodes zi. Within each interval, the Hubble parameter is approximated according to the chosen parametrization, allowing for an analytic expression for ∆dcom,i(z)that accounts for the variation of H(z)across the bin. The total comoving distance to a redshift zcan then be obtained by summing the contributions from all nodes below zand including the partial interval up to zitself: dcom(z) = Nj−1 X j=1 ∆dcom,j(zj+1),(9) where Njis the number of redshift nodes satisfying zj< z. This approach effectively constructs the comoving distance by combining exact contributions within each redshift bin, which is particularly convenient for spectral siren analyses where the source redshifts are sampled discretely. Therefore, Eqs. 4and 9provide all the cosmological ingredients needed to compute the comoving distances required in a spectral siren framework. 3 The spectral siren inference framework 3.1 Hierarchical likelihood and constraint estimators To jointly infer the population and cosmological parameters of our GW catalog, we employ a hierarchical Bayesian inference (HBI) framework [43,44]. In particular, we focus on the spectral 4 IOP Publishing Journal vv (yyyy) aaaaaa Author et al siren approach, which relies solely on GW data, without the use of external galaxy catalogs to estimate source redshifts. The detection of GW signals by terrestrial interferometers can be modeled as an inhomogeneous Poisson process, subject to selection effects arising from the finite sensitivity of the GW detectors. Following [43,44], the hierarchical likelihood of observing Nobs GW events described by a dataset {x}over an observation time Tobs can be written as: L({x}|Λ, H(z)) ∝e−Nexp(Λ) Nobs Y i=1 Tobs ZdθdzL(xi|θ,z,Λ)1 1+z dNcbc(Λ) dθdzdts ,(10) where L(xi|θ, z, Λ)is the likelihood of a single event, with θdenoting the intrinsic parameters of the binary, zthe cosmological redshift, and Λthe population and cosmological hyper-parameters. The correction for selection effects is encoded in the expected number of detections, Nexp(Λ), while the final term in Eq. 10 represents the CBC merger rate. This hierarchical likelihood can also be expressed in an equivalent “scale-free” form, given by: L({x}|Λ, H(z)) ∝ Nobs Y i=1 RdθdzL(xi|θ, z, Λ)1 1+z dNcbc(Λ) dθdzdts Rdθdzpdet(θ, z, Λ, H(z)) 1 1+z dNcbc(Λ) dθdzdts ,(11) where the probability of detection, pdet(θ, z, Λ), accounts for selection effects. For further details on the derivation of this likelihood and the definitions of the various terms, we refer the reader to [14]. In our analysis, the CBC population is described by the set of intrinsic parameters θ= (m1,s, q), where the mass ratio is defined as q=m2/m1. Spins are not considered in this work. With this choice of population parameters, the hierarchical likelihood in Eq. 11 can be rewritten by expressing the CBC merger rate as: L({x}|Λ, H(z)) ∝ Nobs Y i=1 RL(xi|m1,s, q, z, Λ)R(z;Λ)ppop(m1,s, q|Λ)dVc dz 1 1+zdm1,sdqdz Rpdet(m1,s, q, z, Λ, H(z))R(z;Λ)ppop(m1,s, q|Λ)dVc dz 1 1+zdm1,sdqdz ∝ Nobs Y i=1 Zi,(12) where R(z;Λ)denotes the CBC merger rate density as a function of redshift, and dVc dzis the differential comoving volume. In general, the merger rate density can be expressed as the product of the local CBC merger rate, R0, and a redshift-dependent rate function, ψ(z;Λ), which is typically modeled to follow the star formation rate [45]. Importantly, Eq. 12 highlights that the hierarchical likelihood is constructed as the product of the marginal likelihoods of individual events, Zi. From a cosmological perspective, this structure is central: cosmological information enters through the mapping between observed quantities and source-frame properties, encoded in H(z)and the comoving volume element. This allows the spectral siren approach to extract cosmological constraints directly from the observed GW population, without relying on external redshift information. The integrals appearing in the hierarchical likelihood of Eq. 12 are evaluated via Monte Carlo integration [46,47]. This is done using a set of simulated detected GW signals (injections) to estimate the selection function statistically, together with posterior samples from the observed GW events, enabling an efficient numerical approximation of the integrals. To ensure numerical stability and avoid bias in these integrals, we apply two cuts on the effective number of samples: one on the posterior samples and another on the injections, set respectively to 10 and 4Nobs, following the same setup of [28]. Inspired by the method presented in [7], which proposes a way to quantify the contribution of individual GW events to the constraints on the Hubble constant or other population parameters, we construct a similar quantity in this study. Specifically, we define a parameter Ckrepresenting the fraction of GW events contributing to the constraints on the Hubble parameter at a given redshift node H(zk). Accordingly, we define Ckas: Ck=1 Nobs PNobs i=1 |ρ(Zi, H(zk))|Var[Zi] PNobs i=1 Var[Zi],(13) where Zidenotes the marginal likelihood of the i-th event from Eq. 12, and Var[Zi]is its variance. The term ρ(Zi, H(zk)) is the Pearson correlation coefficient between that marginal likelihood and 5 IOP Publishing Journal vv (yyyy) aaaaaa Author et al the Hubble parameter at redshift node zk[48]. This coefficient quantifies the extent to which a particular GW event icontributes to the inference of a given parameter, here H(zk). If this coefficient vanishes for a given event, it implies that variations in the population parameter have no effect on the single-event marginal likelihood. Consequently, the contribution of this event to the population likelihood is constant, and it does not influence the inference of that population parameter. Intuitively, GW events contribute to a given H(zk)if their marginal likelihood correlates with it; the stronger this correlation across the population, the tighter the constraints on the parameter1. In Eq. 13, both the Pearson correlation coefficient and the variance are calculated over the population posteriors. We note that Eq. 13 can be used for any population parameter Λ on which we either set priors or it is implied by the population model, for instance the luminosity distance at a given redshift node. 3.2 Population and cosmological models In addition to the non-parametric model with splines for the Hubble parameter, we must also define the population model for the mass distribution, ppop(m1,s, q|Λ), and the BBH merger rate, R(z;Λ), to construct the hierarchical likelihood in Eq. 12. For the mass distribution, we adopt a model inspired by the fiducial BBH population presented in [49], namely a broken power law + 2 peaks model for the primary mass m1,s, coupled to a smoothed power law for the mass ratio q. The broken power law + 2 peaks model combines a broken power law B, with a slope transition at a break mass mbreak, with two truncated Gaussian peaks G1and G2. The full distribution for the primary mass, π(m1,s|Λ), is then given by p(m1,s|Λ) = (1 −λg)B(m1,s|Λ)+λgλlow gG1(m1,s|Λ)+λg(1 −λlow g)G2(m1,s|Λ),(14) where λgand λlow gare the mixing fractions for the Gaussian peaks. The mass distribution is defined between a minimum mass mmin and a maximum mass mmax. Additionally, two smoothing factors δm,min and δm,max measured in solar masses and representing the width of a sigmoid function are applied to the lower and upper end of the spectrum to encode a possible tapering of the spectrum boundaries. The mass ratio distribution is modeled as a power over the secondary mass and conditioned on the primary mass. Explicitly, the mass ratio distribution π(q|Λ)reads π(q|Λ)∝qβqΘ[mmin/m1,1],(15) For the BBH merger rate, R(z;Λ), describing the evolution of BBH mergers with cosmological redshift, we adopt the commonly used parameterization following the star formation rate from [45]: R(z;Λ) = 1+(1+zp)−γ−κ(1+z)γ 1 + [(1 + z)/(1+zp)]γ+κ,(16) where zprepresents the turning point between two regimes governed by the power indices γand κ. Regarding the cosmological models, we consider two parameterizations: our main non-parametric model, where H(z)is constructed using splines as an autoregressive function in log space (Section. 3), and the parametric flat-ΛCDM model, given by H(z)=H0qΩm,0(1+z)3+(1−Ωm,0),(17) which is commonly used in GW cosmology analyses. For the autoregressive H(z)model, we employ 100 redshift nodes, uniformly distributed in ln(1 + z)up to zmax = 10. The spacing between consecutive nodes i+ 1 and i, denoted ∆ui, is ∆ui= ln(1 + zi+1)−ln(1 + zi)≈0.025,(18) ensuring sufficient resolution in redshift. Between nodes, the Hubble parameter is interpolated using Eq. 4. It follows that evolution of H(z)between two nodes is H(zi+1)=H(zi)1 + 3 2(1+wi∆ui),(19) where wiis the EOS parameter in the i-th redshift bin. We adopt a uniform prior wi∈[−2,1/3], encompassing scenarios from phantom dark energy-dominated universes to radiation-like EOS. 1We note that Ckshould not be interpreted strictly as the fraction of GW events contributing, as this would require a perfectly multivariate Gaussian posterior for Ziand H(zk). However, it provides a useful indication of the most promising redshift ranges in terms of constraining power. 6 IOP Publishing Journal vv (yyyy) aaaaaa Author et al In this non-parametric framework, the cosmological parameters consist of the Hubble parameter at z= 0,ln H0, and the (Nnodes −1) EOS coefficients wi. Table 1summarizes the priors and provides a description for each parameters (including the mass and rate prescriptions). Notably, in the limit ∆ui≪1, the autoregressive model of Eq. 19 approximates the Hubble parameter in Eq. 7. The flat-ΛCDM model is more straightforward in terms of parameterization and is described in detail in [28]. 3.3 Data: GW catalog, injections, and BAO measurements All results presented in this study are based on the recent GWTC-4.0 catalog release by the LVK collaboration, which includes CBCs detected from the first observing run (O1) through the first part of the fourth observing run (O4a) [33]. From the 218 GW candidates, we select a subset of 137 BBH events with a FAR <0.25 yr−1to ensure that the majority of candidates are of astrophysical origin. We also exclude triggers occurring during the engineering run immediately preceding O4a, maintaining consistency with the latest LVK cosmological study [28]. The very massive event GW231123_135430 [50] is also excluded from our BBH sample, as key astrophysical properties, such as the component masses and inferred luminosity distance, appear to be highly sensitive to the waveform model choice during parameter estimation (PE). We emphasize that the event selection in this study is based solely on the FAR threshold. In particular, we use the stable-release version of the GWTC-4.0 catalog, which is also employed in both population and cosmological analyses by the LVK [28,49]. The GW detection sensitivity appearing in the denominator of Eq. 11 is estimated using an injection campaign. Simulated GW signals (injections) are added to the real detector noise and recovered by the LVK search pipelines, and the evolving sensitivity for each observing run is accounted for directly in the injection file. For further details on the production and use of these injections in the Bayesian inference framework, we refer the reader to [51–53]. In addition to the GW-based approach, we complement our analysis by incorporating measurements of the Hubble parameter, H(z), derived from BAOs. We use the latest BAO results from the Dark Energy Spectroscopic Instrument (DESI) [54], based on Data Release 2 (DR2) after three years of observations. BAOs provide a powerful tool for measuring the expansion history, using a characteristic scale imprinted on matter clustering by pressure waves that propagated in the coupled baryon-photon fluid of the pre-recombination Universe. The DESI DR2 sample includes over 14 million galaxy and quasar redshifts, spanning six tracer samples (e.g., luminous red galaxies, emission-line galaxies, quasars). These measurements provide distance constraints in the form of the transverse comoving distance, DM(z)/rd, normalized by the sound horizon at the drag epoch, rd, and the normalized radial distance, DH(z)/rd≡(c/H(z))/rd. For the calculation of H(z), we adopt the value of the sound horizon at the drag epoch as estimated in [1]. 4 Results All the results presented in this article were obtained using Gsirens , a new Python package developed for rapid Bayesian inference with GW data. It relies on Hamiltonian Monte Carlo, specifically the No-U-Turn Sampler (NUTS) implementation in NumPyro and JAX [55,56]. The inference, including the hierarchical Bayesian analysis, was performed on a single NVIDIA GeForce RTX 4050 GPU. With this setup, analyzing the 137 BBH events, using 10,000 samples per event, takes approximately 10 hours with our spline model featuring 100 redshift nodes, and only 2 hours with the commonly used flat-ΛCDM cosmological model. For all the results on the reconstructed mass and mass ratio distributions and BBH merger rate as function of redshifts we obtain consistent results with [28]. In the following, we will focus on the constraints on the Hubble parameter. 4.1 The Hubble parameter with GWs Figure. 2presents the main result of this paper. The left panel shows the maximum a posteriori (MAP) estimation and 68.3% highest density intervals (HDI) of H(z)in each of the i-th redshift bins generated by the autoregressive model. The right panel displays the reconstructed projection of the Hubble parameter, this time assuming the standard parametric flat-ΛCDM cosmological model. In both the parameterized and non-parametric reconstructions, it is clear that GW data encode some level of information about the Hubble parameter, as the inferred values are constrained well within their implied prior ranges. A perfectly non-informative inference would yield broader HDIs, covering the entire prior range, behavior which is not observed here. With both models, we obtain consistent estimates for the local value of the Hubble parameter, in agreement with current 7 IOP Publishing Journal vv (yyyy) aaaaaa Author et al 1 2 3 4 5 1 + z 0 50 100 150 200 250 300 H(z) [km s−1Mpc−1] Prior H(z) Splines 1 2 3 4 5 1 + z Prior H(z) flat-ΛCDM Figure 2. Prior and posterior predictive distributions of the Hubble parameter H(z)as a function of cosmological redshift in log10 space, from a spectral siren inference using 137 GW detections with FAR<0.25 yr−1, from the GWTC-4.0 catalog. Left: Non-parametric inference with the spline H(z)model. Binned constraints are plot as continuous to improve readability. Right: Parametric inference of the Hubble parameter H(z)assuming a flat-ΛCDM cosmological model. The dark colored contours represent the 68.3% highest density intervals (HDI) around the maximum a posteriori (MAP) estimation. literature [1,2,28,35]. We find that the Hubble constant is measured at H0= 85+38 −29 kms−1Mpc−1 and H0= 86+42 −24 kms−1Mpc−1for the non-parametric and parametric approaches, respectively. While the reconstructions of the Hubble parameter are the same at low redshifts, they diverge significantly at higher redshifts. This suggests that the constraints on H(z)and its behavior in this region, are driven primarily by the choice of the assumed cosmological model and its original priors, rather than being data-driven, as one might hope. This effect is further amplified by the fact that beyond a certain redshift (z∼1.2), no GW detections are available, making the shape of H(z)entirely driven by the cosmological model. As a direct consequence, we argue that in spectral siren analyses, it is not trivial to disentangle how much of the observed cosmological constraints truly originate from the GW data itself, as opposed to being a posterior projection induced by the chosen model, in this case, flat-ΛCDM. 4.2 Anchoring H(z)with BAOs To further quantify the induced effects of modeling and posterior projections on the estimation of the Hubble parameter H(z), we extend our analysis by incorporating additional “anchors” that can measure H(z). We include the latest H(z)estimates from BAOs produced by the DESI collaboration, providing five anchor points in the redshift range 0.5<z<1.5. We use the measures of the drag distance DM(z)/rdfrom DESI BAO latest measurements and use the sound horizon measured from the CMB to estimate the Hubble parameter [1]. Hence, we now look at three distinct analyses, all using the non-parametric modeling for H(z): the first relies solely on GW data (identical to the left panel of Fig. 2), the second uses only the five BAO measurements between redshifts 0.5 and 1.5, and the third combines both GW data and BAO anchors within the same inference framework. The results of these inferences are reported in Fig. 3, where the GW measurement is shown in blue, the BAO in purple and the anchored GW with BAO result in yellow. Once again, the MAP is given by the solid line and the colored contours outline the 68.3% HDIs. As shown in Fig. 3, the five BAO data points from DESI provide robust measurements of the Hubble parameter in the redshift range 0.5<z<1.5. Outside this range, the BAO-only reconstruction of H(z)is clearly model-induced, essentially a posterior projection of the priors inherent to our autoregressive model with splines. By anchoring our GW-based H(z)reconstruction with BAO measurements, we first observe no additional constraining power in the redshift range 0.5<z<1.5, as expected given the larger error budget of the GW inference in this region. Below z < 0.5, however, we detect a shift in the Hubble parameter toward lower values, accompanied by a slight reduction in overall uncertainties. This effect arises because this redshift range contains the majority of our GW candidates, while the BAO-only inference here is mostly model-induced, as a projection of the prior. Conversely, above z= 1.5, the combined anchored reconstruction of H(z)aligns again with the BAO-only results indicating that GW data is not informative at z > 1.5. 8 IOP Publishing Journal vv (yyyy) aaaaaa Author et al 1 2 3 4 5 1 + z 50 75 100 125 150 175 200 H(z) [km s−1Mpc−1] GW Splines BAO Splines GW+BAO DESI BAO Figure 3. Reconstruction of the inferred Hubble parameter H(z)using our non-parametric model with splines, for GW detections only (blue), BAOs only (purple), and combining GW and BAO (yellow). The BAO data points estimated by DESI are also shown in dark purple on this figure. For each posterior predictive check shown here, the MAP is represented by the solid colored line and the 68.3% HDI are shown with the contours. 4.3 Quantifying the GW constraining power So far, we have demonstrated two key points: (i) While it is possible to extract cosmological information from GW data alone, beyond a certain redshift, the reconstruction of the Hubble parameter reduces to a prior-induced projection of the population posterior under the assumed cosmological model, whether a parametric or non-parametric model is used. (ii) Anchoring H(z) with external probes such as BAO helps constrain the overall shape of the Hubble parameter and identifies the redshift regions where GWs improve the inferred values, primarily at low redshift where most of the GW candidates lie. Based on the results presented above, one might consider using the non-parametric inference of H(z)at each node to identify the redshift nodes zkwhere GWs provide the tightest constraints on H(z=zk). However, as we have discussed in Sec. 4.2 this can not be easily done by using posterior and prior predictive distributions on H(z), as these are implied by the model. To address this, we now turn to the constraining fraction coefficient introduced in Eq. 13, aiming to pinpoint the redshift bin(s) where GW events are most informative for the Hubble parameter. Here, we emphasize that the optimal redshift bin is not the one providing the tightest constraints on H(z), but the one where the estimate of H(z)is least induced by the model and more data-driven by GWs. Figure. 4shows the estimated constraining fraction coefficients for both the Hubble parameter H(z)(left panel) and the luminosity distance dL(z)(right panel), computed via three approaches: GW-only and GW+BAO with the non-parametric cosmological model; and the strongly parametrized flat-ΛCDM cosmological model. On the left panel, the constraining coefficient differs significantly between the flat-ΛCDM and non-parametric models, particularly in the lower redshift bins. This discrepancy arises because in reality, GW candidates are detected in luminosity distance space, not redshift space where H(z)is defined. Consequently, the constraining power on H(z)also reflects how each model, parametric or non-parametric, and their priors, are mapped from the luminosity distance space. The flat-ΛCDM and non-parametric models thus induce distinct priors on luminosity distance, leading to differing constraints on H(z)across redshift bins. We consequently examine the constraining fraction coefficients in a more natural space, the luminosity distance space, dL(z), as shown in the right panel of Fig. 4. Here, the strong 9