Bright Siren Cosmology Beyond Fisher: DALI Forecasts for Third-Generation Gravitational Wave Observatories

E-mail: [email protected], [email protected], [email protected], [email protected]

Introduction

Gravitational wave (GW) observations have opened an observational window onto the Universe, providing direct access to compact-binary coalescences and enabling distance measurement methods that differ from those obtained through electromagnetic probes. In the standard-siren framework, originally introduced by Schutz [1], compact-binary coalescences can be used as distance indicators because the GW signal directly encodes the luminosity distance $d_L$ to the source, without the need for external calibration [2, 3]. At leading order, the amplitude of the inspiral waveform scales inversely with the $d_L$, but it also depends on the redshifted chirp mass $\mathcal{M}_c$ and on the binary orientation, particularly the inclination angle $\iota$. Therefore, the distance information is not obtained from the amplitude alone.

The time-derivative of the inspiral frequency $\dot{f}$ depends on the $\mathcal{M}_c$ with the same functional dependence as the waveform amplitude. Thus, one concludes that $d_L$ can be determined directly from GW data $h(t)$ and $\dot{f}(t)$. In this sense, compact-binary coalescences, like binaries of black holes (BBH), neutron stars (BNS) or mixed neutron star-black holes (NSBH), behave as absolute distance indicators [1]. In practice, however, the inference of $d_L$ is limited by parameter degeneracies, most notably the $d_L-\iota$ degeneracy, which can generate asymmetric and non-Gaussian posteriors for individual events.

GW observations alone do not directly provide the source redshift $z$, since the waveform depends on detector-frame masses rather than on intrinsic masses and redshift separately. For cosmological applications, an additional redshift measurement is therefore required. In this work, we focus on bright sirens, namely GW events with an identified electromagnetic counterpart. In such cases, the counterpart allows the host galaxy to be identified and its redshift to be measured spectroscopically, so that each event provides a direct point in the $d_L-z$ relationship. Once both $d_L$ and $z$ are available, the source provides a direct and calibration-independent measurement of the cosmic expansion history [2, 3]. In the low-redshift regime, this relation enables an inference of the Hubble constant $H_0$.

The binary neutron-star merger GW170817 provided the first GW standard siren measurement of the Hubble constant [4]. Electromagnetic observations identified an optical transient spatially associated with the galaxy NGC 4993, allowing the host-galaxy redshift to be measured independently of the GW signal. By combining $d_L$ inferred from the GW data with the recession velocity inferred from electromagnetic observations, it was obtained $H_0 = 70.0^{+12.0}_{-8.0}\ \mathrm{km\ s^{-1}\ Mpc^{-1}}$ [4].

It has been shown that as the number of detections increases, our constraints on $H_0$ can become more precise, see [5, 6, 7]. Chen et al. [5], for instance, evaluated the ability of second-generation GW detectors to constrain the $H_0$ using both BNS and BBH mergers. Redshift estimation was based either on electromagnetic counterparts, as in GW170817, or on the redshift distribution of host galaxies or galaxy clusters. Their results indicate that observations of approximately 10, 60, and 200 BNS mergers with electromagnetic counterparts could constrain $H_0$ to precisions of 4%, 2%, and 1%, respectively. Unfortunately, this prediction is limited by the small number of observed BNS mergers [8].

Third-generation GW observatories (3G), like Einstein Telescope (ET) and Cosmic Explorer (CE), will provide a substantial advance over current second-generation detectors, particularly in two aspects: the sky localization of nearby compact-binary sources and the precision for which $d_L$ can be inferred [9, 10, 11]. Better localization increases the likelihood of associating GW events with electromagnetic counterparts, while more accurate distance measurements reduce one of the dominant uncertainties in standard siren estimates of $H_0$. Zhao et al. [12] investigated different configurations of 3G detector networks and evaluated their ability to constrain cosmological parameters using BNS and NSBH mergers. They found that a network of ET-like detectors distributed across the United States, Europe, and Australia could measure the $d_L$ of a BNS merger at $z = 0.1$ with a fractional uncertainty ranging from about 10% to 1%. For a sample of 100 comparable events with identified electromagnetic counterparts at that redshift, the resulting measurement of the $H_0$ could reach a precision of approximately 0.3%.

This expectation is supported by several independent forecasts. Taylor and Gair [13] found that ET could catalog tens to hundreds of thousands of BNSs per year, and that a sample of $10^5$ detections could constrain dark-energy equation-of-state parameters with an accuracy comparable to future CMB+BAO+SNIa measurements. Similarly, BBHs forecasts for 3G networks indicate that an ET+CE+CE configuration, combined with a galaxy catalog information complete up to $z = 1$, can constrain $H_0$ to about 0.8% and $\Omega_m$ to about 10% within one year of observations. On the other hand, bright siren forecasts with future ET and CE observations further indicate that BNS and NSBH events with electromagnetic counterparts can provide independent constraints on the redshift evolution of the dark-energy equation of state [14, 15]. Together, these forecasts indicate that with 3G detectors, GW astronomy will become a precise probe of the late-time expansion of the Universe.

However, robust forecasting in this context requires statistical methods that remain reliable when the likelihood departs from Gaussianity. For this reason, Fisher-matrix (FM) methods, although computationally efficient and widely used [16, 17, 18, 19, 20], are not always accurate when the likelihood exhibits strong curvature, skewness, or non-linear correlations [21, 22, 23, 24, 25, 26]. Higher-order likelihood expansions provide an extension of the Fisher approach while retaining its derivative-based structure. In this work, we investigate this problem using the Derivative Approximation for Likelihoods (DALI) formalism [21, 22, 23, 27]. We compare the Gaussian Fisher approximation with higher-order DALI expansions in their ability to re-construct the likelihood of bright standard sirens and to propagate the inferred distance information to cosmological parameter estimation.

Although standard sirens can also be used statistically in the absence of an electromagnetic counterpart, such dark-siren approaches are not the focus of the present analysis. Here we restrict attention to sources with electromagnetic counterparts, for which the redshift information is directly available, and the central statistical challenge relies on understanding how non-Gaussian structures in the GW likelihood propagate into cosmological-parameter inference. Therefore, bright sirens provide the setting adopted in this work to assess the performance of higher-order likelihood expansions beyond the Fisher approximation.

Our analysis is based on simulated bright siren catalogs for 3G networks. We assess how accurately Fisher and DALI approximations reproduce the underlying likelihood structure, with particular attention to the $d_L - \iota$ degeneracy, and we examine the consequences of these approximations for cosmological inference. Our goal is to determine under which conditions the Gaussian approximation remains adequate and when higher-order DALI corrections become necessary for reliable forecasting with future GW observations.

Recent studies of the ET science case have investigated the cosmological potential of bright standard sirens under different assumptions about detector configurations and electromagnetic follow-up strategies. For a 10-km triangular ET configuration, forecasts report marginalized uncertainties on the Hubble constant of approximately 0.9% for a one-year gravitational-wave catalog with kilonova counterparts and 5.7% for a five-year catalog with gamma-ray burst counterparts [28, 29]. These results illustrate how the expected cosmological precision depends on the assumed source population, detector configuration, and electromagnetic follow-up strategy. In this context, rather than providing an updated end-to-end forecast for the ET science case, the present work investigates how non-Gaussian gravitational-wave likelihood approximations affect cosmological inference for a fixed bright-siren catalog, focusing on the accuracy of Fisher and DALI-based methods.

The remainder of this paper is organized as follows. Section 2 reviews the Fisher matrix and DALI formalisms, highlighting their underlying assumptions and limitations. Section 3 describes the construction of the simulated gravitational-wave catalogs, while Section 4 presents the gravitational-wave parameter-estimation methodology, including the Fisher and DALI approximations. Section 5 presents the resulting cosmological constraints and provides a quantitative comparison of the different approximations. Finally, Section 6 summarizes our main findings and discusses prospects for future work.

Likelihood Approximations

Throughout this work, we refer to the exact likelihood as the unapproximated functional form used for the inference, namely the standard Gaussian-noise likelihood in the GW analysis and the $\chi^2$-based likelihood in the cosmological analysis. We also work with geometrized units, where $G=c=1$.

Fisher matrix formalism

The Fisher information matrix provides a convenient way to predict the expected uncertainties on model parameters before any data are collected. For a model with $\mathcal{N}$ parameters, it typically requires only $O(\mathcal{N})$ likelihood or model evaluations [23]. This scaling should be understood as follows: the Fisher matrix is constructed from derivatives of the theoretical predictions with respect to each parameter, evaluated around a fiducial point [30, 31]. In a finite-difference implementation, one evaluates the model at the fiducial point and at a fixed number of displaced points along each of the $\mathcal{N}$ parameter directions. For instance, a central-difference scheme requires approximately $1+2\mathcal{N}$ evaluations. Since constant factors are ignored in asymptotic notation, this cost scales linearly with the number of parameters, i.e. as $O(\mathcal{N})$. Once these derivative vectors are obtained, the individual Fisher-matrix elements are computed through algebraic contractions, without requiring a full exploration of the likelihood surface.

The method follows from a Taylor expansion of the log-likelihood around its maximum. The derivatives of the theoretical predictions with respect to the parameters quantify how sensitive the observables are to small changes in the parameter vector $\boldsymbol{\theta}$. If the likelihood $\mathcal{L}(\boldsymbol{\theta})$ is well approximated by a Gaussian in the neighborhood of its peak, then $\mathcal{L}$ is locally quadratic and the parameter constraints are determined by the curvature of $\log\mathcal{L}$ [30, 31], which can be written as

where $\mathbf{F}$ is the Fisher information matrix and $\Delta\boldsymbol{\theta}$ represents deviations from the best-fit parameter values. In this approximation, the inverse of the Fisher matrix corresponds to the covariance matrix of the estimated parameters,

such that

gives the minimum expected uncertainty for each parameter. This result is known as the Cramér–Rao bound.

In the GW data analysis context, codes such as GWBENCH [17], GWFAST [20] and GWFISH [18] implement this approach to rapidly evaluate the sensitivity of proposed observatories (e.g., ET, CE, LISA) or multi-band networks. Similarly, in cosmology, the Fisher information matrix is widely used to forecast constraints on parameters such as the matter density parameter $\Omega_M$, the amplitude of matter fluctuations $\sigma_8$, the dark-energy equation-of-state parameter $w_0$, or the sum of neutrino masses [32].

Nevertheless, Fisher forecasting has important limitations. The Fisher matrix method provides an exact description only when the posterior is Gaussian in the parameters; consequently, it predicts symmetric, elliptical constraint contours. In practice, this approximation can break down, particularly in the presence of strong parameter degeneracies or when the data are only weakly informative, leading to inaccurate uncertainty estimates [23].

Derivative Approximation for Likelihoods (DALI)

If the posterior is non-Gaussian, the Fisher approximation can become inaccurate, particularly when the observables depend nonlinearly on the model parameters, when the likelihood displays curved degeneracy directions, or when weakly constraining data produce non-elliptical posterior contours [23, 24, 26, 33, 34]. To overcome these limitations, the Derivative Approximation for Likelihoods (DALI) extends the Fisher formalism by incorporating higher-order derivatives of the likelihood or, equivalently, of the theoretical model, while preserving a positive definite and normalizable approximation to the posterior [22, 23, 27, 34, 35]. In its lowest non-Gaussian correction, known as the Doublet approximation, the method requires only $\mathcal{O}(\mathcal{N}^{2})$ likelihood or model evaluations for $\mathcal{N}$ parameters [23, 27].

This scaling arises because the Doublet approximation depends not only on the first derivatives of the theoretical observables, $\boldsymbol{\mu}_{,a}$, but also on their second derivatives, $\boldsymbol{\mu}_{,ab}$. These second derivatives must be evaluated for all independent pairs of parameters $(a,b)$. For a model with $\mathcal{N}$ parameters, the number of such independent pairs is

which scales quadratically with the number of parameters. In a finite-difference implementation, each derivative is obtained from a fixed number of model evaluations around the best-fit point. Hence, although the numerical prefactor depends on the derivative stencil adopted, the leading computational cost scales as

Thus, Doublet is more expensive than the Fisher approximation, whose leading cost is associated with $\mathcal{O}(\mathcal{N})$ first-derivative evaluations, but it remains much less costly than a full sampling of the exact likelihood [23, 27, 36].

Following the same local-expansion logic underlying Fisher forecasting, let us define

where $\hat{\theta}^{a}$ denotes the best-fit. Since the first derivatives vanish at the maximum, the Fisher approximation corresponds to retaining only the quadratic term of the expansion,

with

On the other hand, DALI reorganizes the expansion not only by powers of $\Delta\theta^{a}$, but also by the order of derivatives of the model predictions $\boldsymbol{\mu}$ [35]. For a likelihood given by:

where $d$ is the vector of data and $M$ encodes the uncertainties of the measurements, the second-order derivative expansion, known as Doublet, is

where commas denote partial derivatives and multiple upper indices in $\Delta\theta$ denotes multiples products, e.g. $\Delta\theta^{ab}\equiv\Delta\theta^{a}\Delta\theta^{b}$, $\Delta\theta^{abc}\equiv\Delta\theta^{a}\Delta\theta^{b}\Delta\theta^{c}$ and so on.

The Triplet approximation extends this construction by including third derivatives of the model predictions:

Thus, Fisher contains only the first derivatives of the model and produces a Gaussian approximation in parameter space; Doublet includes second derivatives, $\mu_{,ab}$; and Triplet includes third derivatives, $\mu_{,abc}$. The method reduces to the Fisher limit when the posterior is effectively Gaussian, but it can capture curved, skewed, or otherwise non-elliptical posterior contours when the model depends nonlinearly on its parameters (see Figure 1) [23, 27, 34, 35, 36].

In GW data analysis, the data vector is replaced by the measured strain $g(t)$, and the theoretical prediction is the waveform $h(t;\theta)$. Assuming stationary Gaussian noise with one-sided power spectral density $S_n(f)$, the likelihood can be written as

where the noise-weighted inner product is defined as

where $N_{\mathrm{dets}}$ is the number of GW detectors.

With the replacement

the Fisher matrix approach becomes

For the Doublet approximation we have [34]

Finally, the Triplet approximation for GW becomes [34]

Here we work with the zero-noise approximation of GW data analysis, where the best-fit parameters $\hat{\theta}$ will be just the fiducial parameters used in the simulations.

Therefore, in the GW formulation, the Fisher approximation is built from products of first waveform derivatives, Doublet adds the terms involving second waveform derivatives, and Triplet further incorporates third waveform derivatives. This structure is particularly useful when the Fisher matrix becomes nearly singular or non-invertible, since higher-order derivative terms can provide a more faithful local approximation to non-Gaussian parameter degeneracies than the purely Gaussian Fisher approximation. In this sense, DALI offers an intermediate approach between analytic Fisher forecasts and full sampling of the exact likelihood [21, 23, 34]. The GW implementation of this formalism was developed in the public Python package GWDALI [21], which computes Fisher, Doublet, and Triplet approximations for arbitrary networks of terrestrial GW detectors and may interface with waveform models from LALSuite and stochastic sampling through Bilby [21]. Another software which make use of this formalism is the hybrid Mathematica/Python-based code called SymDALIwhich make use of symbolic derivatives to compute DALI tensors performing stochastic sampling via Bilby [37]. Ryan et al. [23] adapted the same DALI philosophy to cosmological forecasting, using the difference between Fisher and DALI likelihood contours as a diagnostic for the validity of Fisher forecasts. Other developments have extended DALI to parameter-dependent covariance matrices and to variational posterior reconstructions for expensive likelihoods [35, 36].

In the next section we simulate our GW catalogs in order to build Hubble-diagrams for cosmological inferences.

Building the data

The mock catalogs employed in this analysis were constructed using GWDALI v0.1.4. The software was used to generate the simulated source populations, calculate the detector responses, impose the detection criterion, and evaluate the Fisher and DALI expansions considered in the subsequent analysis. The physical assumptions entering the simulations, including the source-population distributions, waveform model, detector network, and fiducial cosmology, are specified independently below.

We generate 1000 statistically independent mock catalogs, each containing 300 detected GW events with electromagnetic counterparts with network signal-to-noise ratio satisfying $\rho_{\mathrm{net}} > 8$. Because an electromagnetic counterpart is assumed to be identified for every selected event, the source redshift is treated as directly measured. These catalogs therefore represent bright standard sirens and provide a controlled setting in which the injected source parameters are known. This allows the accuracy and limitations of the Fisher and DALI likelihood approximations to be assessed directly.

Source redshift distribution.

The sources are distributed over the redshift interval

The adopted distribution is intended to represent the observable population of bright sirens expected from 3G detectors operating jointly with a dedicated electromagnetic mission such as THESEUS. Considering the parametrization adopted in [38, 39], the redshifts of the sources are sampled from the following normalized distribution

where $N_z$ is the normalization constant.

Binary population and angular distributions.

The simulated sources are equal-mass compact binaries satisfying

The detector-frame masses are obtained from the corresponding source-frame masses through

The sky coordinates of the sources ($\alpha$ and $\delta$) are isotropically sampled over the celestial sphere. The polarization angle is sampled uniformly over its appropriate angular interval. These choices ensure that the simulated population does not contain a preferred sky direction.

To approximate the selection effects associated with the observation of an electromagnetic counterpart, the binary inclination angle, $\iota$, is drawn from a smooth-cutoff distribution. Following [40], we adopt

where

Throughout this work, we adopt

The parameter $\iota_0$ acts as an effective angular scale for the electromagnetic selection of bright sirens. Joint GW–GRB observations preferentially select binaries viewed close to the orbital angular-momentum axis, corresponding to nearly face-on configurations ($\iota\simeq0$) or nearly face-off configurations ($\iota\simeq\pi$) [41, 42]. In this regime, the GRB viewing angle is closely related to the binary inclination, and the identification of a GRB counterpart therefore favors small values of $\iota$ or $\pi-\iota$.

This selection is particularly relevant because the luminosity distance inferred from the GW amplitude is strongly correlated with the inclination angle [42, 43]. Restricting the population to nearly face-on or face-off systems partially reduces this degeneracy and improves the expected luminosity-distance measurements.

The value $\iota_0=20^\circ$ is consistent with previous standard-siren forecasts in which angular scales of a few tens of degrees are used to model GRB beaming and electromagnetic selection. Nissanke et al. [43], for example, considered a beamed short-GRB population in which most detectable systems were observed within an angle of approximately $25^\circ$. Similarly, Belgacem et al. [44] adopted a GRB beaming angle of approximately $25^\circ$ in mock catalogs of joint GW–GRB standard sirens. The adopted value therefore provides a smooth analytic approximation to the angular preference for bright sirens viewed close to the binary angular-momentum axis, rather than a complete model of electromagnetic detection.

The smooth-cutoff inclination distribution adopted here provides an analytic approximation to the angular selection associated with electromagnetic-counterpart detection, avoiding the abrupt transition implied by a hard angular cutoff. However, this prescription does not explicitly model the GRB emission, the instrumental selection function of a specific satellite, or the probability of obtaining a redshift for each event. Instead, we assume that an electromagnetic counterpart and its redshift are available for every source included in the simulated catalogs. A more detailed treatment of these observational effects would be a valuable extension of this work.

Figure 2 summarizes the source-population model used in the simulations. The left panel shows the adopted redshift distribution, while the right panel illustrates the inclination selection associated with the assumed electromagnetic-counterpart population.

Gravitational Wave signal model

The GW signals are modeled in the frequency domain using the TaylorF2 waveform approximant. TaylorF2 is obtained through the stationary-phase approximation and describes the inspiral phase of a quasicircular compact binary. In the restricted post-Newtonian approximation, the amplitude is retained at leading order, while the Fourier phase is included up to 3.5 post-Newtonian order [45].

The two gravitational-wave polarizations in the Fourier domain are given by

where $v \equiv(\pi\mathcal{M}_z f)^{1/3}$ is the characteristic orbital velocity of the binary, $\mathcal{M}_z$ denotes the detector-frame total mass, and $t_c$ and $\phi_c$ are the coalescence time and phase, respectively.

The post-Newtonian (PN) coefficients $\varphi_k$ can be found in [46]. Here, $d_L$ is the luminosity distance, $\iota$ is the inclination angle of the binary, and $\mathcal{M}_{c,z}$ is the detector-frame chirp mass, related to the source-frame chirp mass by

where $M$ is the source-frame total mass and $\eta$ is the symmetric mass ratio.

Detector network and detection criterion.

The simulated detector network consists of ET and CE. CE is modeled as a single L-shaped interferometer, with an opening angle of $90^\circ$ and the sensitivity curve provided in [47]. ET is modeled in its triangular configuration, equivalent to three co-located interferometers with opening angles of $60^\circ$, rotated by $0^\circ$, $120^\circ$, and $-120^\circ$, respectively. For ET, we adopt the historical ET-D sensitivity curve available in [48], as implemented in GWDALI¹. The network therefore comprises four detector channels in total.

Detector channel	Observatory	Longitude	Latitude	Rotation	Arm opening angle
det0	CE	$-119^\circ$	$46^\circ$	$45^\circ$	$90^\circ$
det1	ET	$10^\circ$	$43^\circ$	$0^\circ$	$60^\circ$
det2	ET	$10^\circ$	$43^\circ$	$120^\circ$	$60^\circ$
det3	ET	$10^\circ$	$43^\circ$	$-120^\circ$	$60^\circ$

For detector “$a$”, the optimal matched-filter signal-to-noise ratio is defined as

where $S_{n,a}(f)$ is the one-sided noise power spectral density of detector $a$. Here, $\tilde{h}_a(f)$ denotes the GW signal measured at the output of detector $a$, obtained as a linear combination of the two GW polarizations:

where $F_{+,a}$ and $F_{\times,a}$ are the antenna-pattern functions, which describe the directional response of detector $a$ to the two GW polarizations. They are given by [46]:

Here, $\Omega_a$ denotes the opening angle between the interferometer arms, while $\alpha_a$ and $\beta_a$ specify the source sky position in spherical coordinates defined with respect to detector $a$. The angle $\psi_a$ is the GW polarization angle in the detector’s reference frame. Finally, $\tau_a$ denotes the GW arrival-time delay at detector $a$ relative to a chosen reference frame, such as the Earth’s center or another detector.

Thus, the network signal-to-noise ratio is

Candidate sources are generated until each catalog contains 300 events satisfying the detection condition $\rho_{\mathrm{net} > 8}$.

Fiducial cosmology

The luminosity distance associated with each simulated redshift is calculated assuming a spatially flat $\omega_0\omega_a$CDM cosmology. The dark-energy equation of state is parametrized as

For a flat Universe, the dimensionless expansion rate is

The fiducial cosmological parameters used to generate the catalogs are

GW Inference

The GW analysis applies the GWDALI package [21]. Following the specifications introduced in the previous section, the CBC signals are modeled in the frequency domain using the TaylorF2 approximant. The set of free GW parameters entering the analysis is

where $\psi$ is the polarization angle, and $\phi_{\mathrm{coal}}$ is the coalescence phase. The quantity $X$ denotes the luminosity-distance variable and is taken to be either $d_L$ or $1/d_L$, depending on the adopted inference prescription.

Since the Fisher matrix provides a local quadratic approximation to the likelihood, the way in which this information is converted into a luminosity-distance measurement can have a non-negligible impact on the resulting GW catalogue. We therefore consider several Fisher-based prescriptions, differing in the distance parametrization, matrix inversion, inclusion of prior information, and posterior construction.

Fisher-Based Approaches

We consider three different prescriptions for the Fisher information. The first corresponds to the original GWDALI implementation used as our baseline, while the remaining approaches were introduced to investigate the dependence of the inferred luminosity-distance uncertainties on the treatment of poorly conditioned Fisher matrices and on the adopted prior information.

Direct Fisher inversion. We construct the Fisher matrix using the inverse luminosity distance, $y \equiv1/d_L$, and subsequently transform it to the $d_L$ parametrization before inversion. Since

the Fisher elements involving the distance parameter transform as

where $\theta$ denotes any of the remaining GW parameters. The resulting matrix is rescaled by its diagonal elements to improve numerical conditioning and is then inverted using the ordinary matrix inverse,

The reliability of each inversion is explicitly assessed through

and only events satisfying $\epsilon_{\mathrm{inv}} < 0.05$ are retained. Singular matrices, matrices with invalid diagonal elements, and covariance matrices with non-finite or negative diagonal elements are also rejected.

Fisher matrix with Gaussian priors. In this prescription, the same transformation to the $d_L$ parametrization is first applied. Independent Gaussian priors are then introduced to regularize poorly constrained directions in parameter space. Their contribution to the Fisher matrix is written as

with prior widths

following the parameter ordering $\{d_L,\iota,\psi,\phi_{\rm coal},\mathcal{M}_c,\eta\}$. The total information matrix is therefore

The covariance matrix is obtained using the same diagonal rescaling and direct-inversion procedure adopted above,

and the same criterion $\epsilon_{\rm inv}<0.05$ is applied to assess the numerical reliability of the inversion. For the accepted events, the Fisher posterior is represented by a multivariate Gaussian centered on the injected source parameters,

from which $10^4$ samples are directly generated. No MCMC sampling is employed in this procedure.

Singlet-DALI. Following [38], a conceptually different prescription is adopted in which the Fisher matrix is used to define a Gaussian approximation to the event likelihood,

where $\Delta\boldsymbol{\theta}=\boldsymbol{\theta}-\boldsymbol{\theta}_0$ and $\boldsymbol{\theta}_0$ denotes the injected parameters. The corresponding posterior is then obtained through

and explored numerically using Monte Carlo sampling. In contrast to the previous approaches, the posterior samples are obtained by sampling the Fisher likelihood together with the adopted priors, rather than being drawn directly from the covariance matrix. This procedure is sometimes referred to as Fisher sampling, MCMC Fisher, or Singlet-DALI [50].

Doublet-DALI

Beyond the Fisher approximation, we also consider the Doublet-DALI expansion. For the Doublet calculation, we adopt $1/d_L$ as the distance variable. This reparametrization is motivated by the behavior of the approximate likelihood in the $(d_L,\iota)$ sector: when the expansion is performed directly in $d_L$, the Doublet reconstruction can exhibit a pronounced bimodality associated with the luminosity-distance–inclination degeneracy. Expressing the likelihood in terms of $u=d_L^{-1}$ leads to a smoother local structure and improves the numerical stability of the approximation [50], as illustrated in Figure 4.

A summary of the priors adopted in the Doublet GW analyses is presented in Table 2. The priors for $\mathcal{M}_c$, $\eta$, $\psi$, $\phi_{\mathrm{coal}}$, and $\iota$ are common to both analyses, while the distance-related prior is implemented directly in $d_L$ for the Fisher analysis and in $u \equiv d_L^{-1}$ for the Doublet analysis. For the inclination angle, we adopt a conservative truncated prior, motivated by the fact that the population of bright sirens is expected to be affected by the selection of events with detectable electromagnetic counterparts and therefore need not follow the isotropic inclination distribution of the underlying binary population. Rather than adopting the inclination distribution of the simulated injections as the prior, we use the same broad, truncated prior for all events. This choice avoids encoding the specific inclination distribution assumed in the simulations directly into the inference and reflects the fact that the inclination distribution of the observed bright-siren population is not known a priori.

Parameter	Prior range	Prior type
$\mathcal{M}_c$	$\mathcal{M}_{c,*}(1-10^{-3}) \leq\mathcal{M}_c \leq\mathcal{M}_{c,*}(1+10^{-3})$	Uniform
$\eta$	$\eta_*(1-10^{-3}) \leq\eta\leq\eta_*(1+10^{-3})$	Uniform
$\psi$	$[0,\pi]$	Uniform
$\phi_{\mathrm{coal}}$	$[-\pi,\pi]$	Uniform
$\iota$	$[0,\pi]$	Sine with cutoff between $[45^\circ,135^\circ]$
$d_L$	$[4\cdot10^{-3},46.7]$ Gpc $\sim z \in[0.01,5.0]$	Comoving-volume prior
$u \equiv d_L^{-1}$	$u \in[1/d_{L,\max},1/d_{L,\min}]$	Transformed comoving-volume prior

Posterior samples are then obtained with a nested-sampling algorithm and are used to characterize the uncertainty associated with the distance-related parameter inferred from each event. The main numerical inputs used in GWDALI are listed in Table 3.

Input	Value used
approximant	TaylorF2_py
sampler_method	nestle
npoints	300
step_size	$10^{-5}$
rcond	0
detectors	det0, det1, det2, det3

Hubble Diagram from GWs via Fisher/DALI

Fisher Enhancements for GW

The different Fisher-based prescriptions described in § 4.1 were first compared at the level of the reconstructed Hubble diagram. Figure 5 shows the luminosity-distance measurements obtained with the direct Fisher inversion, the Fisher matrix with Gaussian

priors, and the Singlet approach. The measurements and uncertainties shown in the left panel are averaged over the 1000 simulated catalogs, allowing the systematic behavior of each prescription to be compared over the full redshift interval.

The direct inversion of the Fisher matrix produces the largest luminosity-distance uncertainties among the considered prescriptions. Although only events satisfying the inversion criterion $\epsilon_{\mathrm{inv}} < 0.05$ are retained, the marginalized uncertainties increase rapidly with redshift. This behavior follows from the poorly constrained directions of the GW parameter space that remain present after the rejection of numerically unreliable inversions, particularly those associated with the luminosity-distance–inclination correlation.

Adding Gaussian prior information to the Fisher matrix modifies the treatment of weakly constrained parameter directions before inversion. The resulting uncertainties are smaller than those obtained with the direct Fisher inversion. The Gaussian priors therefore regularize the covariance estimation without reproducing the golden-event population observed in the pseudo-inverse prescription.

The Singlet approach differs from the previous prescriptions because the Fisher matrix is not inverted and interpreted directly as the final posterior covariance. Instead, the quadratic Fisher likelihood is sampled together with the adopted physical priors. Among the Fisher-based treatments considered here, this prescription produces the luminosity-distance reconstruction closest to that obtained from the exact GW likelihood. The right panel of Figure 5 compares the Hubble diagram obtained from the ensemble-averaged Singlet catalogs with the exact reconstruction from a representative catalog, showing consistent distance measurements over the analyzed redshift range.

These results demonstrate that the construction of a Fisher-based Hubble diagram depends not only on the quadratic approximation to the GW likelihood, but also on how this approximation is converted into marginalized parameter constraints. All approaches lead to distinct luminosity-distance uncertainties even when they originate from the same local Fisher information. Among these treatments, Fisher sampling provides the closest agreement with the exact GW reconstruction and is consequently adopted as the reference Fisher enhancement in the subsequent comparisons.

The comparison with a single GW-Exact realization shown in Figure 5 is intended only as an illustrative reference, providing a visual indication of how the mean reconstructed quantities from the approximate catalogs compare with an individual exact realization. A systematic validation of the GW Fisher and DALI approximations against the exact GW likelihood was presented in our previous work on the validation and release of GWDALI v1.0 [50]. Here, we focus instead on the propagation of the GW approximations to cosmological inference across the simulated catalog ensembles.

Doublet Approach

The reconstructed mean distance–redshift relation shown in Figure 6 reveals a qualitatively different behavior between the Fisher (different approaches) and Doublet approximations. The Fisher reconstruction remains approximately centered on the fiducial relation across the full redshift range, although its uncertainty band broadens rapidly with increasing redshift. In contrast, the Doublet reconstruction produces substantially narrower uncertainty bands, but its mean is systematically displaced toward smaller luminosity distances.

This comparison therefore highlights a trade-off between precision and accuracy in the present implementation. While the Doublet approximation better constrains the local shape of the likelihood, the resulting reduction in posterior width is accompanied by a bias in the recovered mean distance. Since the DALI tensors are computed through automatic differentiation, this behavior is more plausibly associated with the truncation of the expansion at the Doublet level, for which a local higher-order approximation may still be insufficient to reproduce the global structure of a strongly non-Gaussian posterior.

The discrepancy becomes more pronounced at higher redshift, where the inference progressively departs from the high signal-to-noise regime in which the Fisher approximation is expected to perform best [34]. In this regime, the $d_L$–$\iota$ degeneracy becomes increasingly important, making the inferred luminosity distance more sensitive to differences in the local likelihood approximation. The recovery of $d_L$ also depends on the detector geometry, the network capability to resolve the GW polarizations, and the adopted prior on $\iota$ [40]. The downward displacement observed for the Doublet reconstruction can therefore be interpreted as a consequence of applying a truncated local expansion to an inference problem with significant non-linear parameter correlations.

Taken together, these Figure 5, Figure 6 indicate complementary limitations of the two approximations. The Fisher matrix fails to adequately capture the joint $(d_L,\iota)$ structure in the present scenario, as its purely quadratic form leads to artificial angular singularities and increasingly large uncertainties at high redshift. The Doublet approximation, in turn, provides a more regular description of the likelihood geometry and significantly reduces the inferred uncertainties, but introduces a redshift-dependent bias in the recovered luminosity distance. In this context, DALI can be regarded as a controlled extension of the Fisher formalism, capable of capturing relevant non-Gaussian corrections and elucidating the limitations of a purely ellipsoidal approximation to the likelihood [23, 27].

Cosmological Constraints

Cosmological Model

The cosmological analysis treats the simulated GW detections as bright standard sirens. The redshift of each event is assumed to be known exactly from the electromagnetic counterpart and is taken directly from the mock catalogs. The luminosity-distance information and its associated uncertainty are inherited from the GW analysis described above. The cosmological parameter vector is defined as

corresponding to the CPL parametrization of the dark-energy equation of state ($\omega_{DE}(a) = \omega_0 + (1-a)\omega_a$ where $a$ is the scale factor). The Hubble expansion rate is written as

and the luminosity distance is given by

For a catalog of GW sources, the cosmological log-likelihood is taken to be

where $d_{L,n}^{\mathrm{obs}}$ and $\sigma_{d_L,n}$ denote, respectively, the luminosity distance inferred from the GW reconstruction and its associated uncertainty.

Semi-analytic Inference

To understand how the measurements of $d_L$ directly affect the measurement of $H_0$, we approximate the luminosity-distance relation by its first-order cosmographic expansion in redshift, as follows:

Let us also to approximate our likelihood via Fisher approximation. In this case, we have a single 1D Fisher Matrix given by:

leading to the uncertainty in $H_0$ as:

Hence, when we have a single GW event, the uncertainty in $H_0$ has a beautiful relation with respect to $d_L$ given by:

Such a result implies tell us that the better the distance measurements, the better the $H_0$ measurements.

Numerical Inference

For the cosmology case, the likelihood is approximated through derivative expansions around the best-fit, calculated. The derivatives of the luminosity-distance function with respect to the cosmological parameters are evaluated using automatic differentiation (using the Python package JAX [51]), rather than finite-difference schemes, ensuring machine-precision derivatives and avoiding the need to choose numerical step sizes.

The approximate cosmological likelihood is sampled with Bilby, adopting uniform priors on all cosmological parameters. The priors are taken to be broad enough, so that they do not artificially truncate the posterior. The posterior samples were obtained with the nestle nested-sampling algorithm implemented in Bilby, using npoints=1000. The uniform prior ranges adopted for the cosmological parameters are listed in Table 4, where $H_0$ is given in km s$^{-1}$ Mpc$^{-1}$.

Parameter	Prior	Range
$H_0$	Uniform	$[50,90]$
$\Omega_m$	Uniform	$[0.01,0.9]$
$\omega_0$	Uniform	$[-5,5]$
$\omega_a$	Uniform	$[-10,10]$

The cosmological inference is performed cumulatively as a function of redshift. For this purpose, the total interval $z \in[0,1]$ is divided into ten cumulative redshift intervals, with upper limits $z_{\max}=0.1,0.2,\ldots,1.0$. For each maximum redshift $z_{\max}$, we select all sources with $z_i < z_{\max}$, so that the number of events increases as the redshift upper limit is raised. The cumulative catalogs used at a given redshift threshold is therefore defined as $\mathcal{D}(z_{\max}) = \left\{d_{L,i}^{\mathrm{obs}}, \sigma_{d_L,i}, z_i \mid z_i < z_{\max}\right\}$. Thus, the inference at $z_{\max} = 0.1$ uses only the nearest sources, while higher redshift thresholds include all previously selected events plus the additional sources located between the previous and current redshift cuts.

In addition to the Fisher and Doublet approximations, we also perform the cosmological inference using the exact likelihood for the luminosity distance data, following equation (55). The Fisher and Doublet likelihoods are instead obtained from the corresponding local derivative expansion around the peak of marginalized posteriors. For each value of $z_{\max}$, the exact, Fisher, and Doublet likelihoods are sampled independently with Bilby using “nestle” with nlive=300.

Uncertainties, Biases and Divergences

Along this work we refer the bias as the displacement of some parameter measurement $\langle X\rangle$ from its fiducial value $X_{fid}$ weighted by the uncertainty $\sigma_X$ on that measurement (with a confidence level of 68.27%)

Following [50] we quantify the comparisons among posteriors (exact against Fisher/DALI approximations) by using the statistical measurement called Jessen-Shannon Divergence (JSD) which quantifies how far some distribution $Q$ is from another distribution $P$ where mathematically is given by [52]:

which is defined in terms of the Kullback–Leibler divergence (KLD) defined below [53].

Cosmological Results

We now compare the propagation of the three GW inference prescriptions, Fisher Singlet, Fisher with Gaussian Priors and Doublet, into the cosmological analysis. For each catalog construction, the cosmological likelihood is independently evaluated using the exact, Doublet, Singlet, and Fisher approximations. This distinction separates the approximation used to reconstruct the GW luminosity-distance information from the approximation subsequently used to infer the cosmological parameters.

The comparison therefore addresses two related effects: how the GW-level likelihood approximation modifies the information entering the Hubble diagram, and how the cosmological likelihood approximation propagates this information into constraints on $H_0$, $\Omega_m$, $\omega_0$, and $\omega_a$. Across the three GW catalog constructions, the cosmological Doublet approximation follows the exact inference more closely than the Singlet and Fisher approximations for most parameter–redshift combinations, as also quantified by the Jensen–Shannon divergences discussed above.

Figure 8–Figure 10 show that the cosmological constraints depend on both the GW catalog construction and the likelihood approximation adopted at the cosmological stage. The differences among the GW-Fisher Gaussian, GW-Singlet, and GW-Doublet catalogs are propagated into the recovered central values and credible intervals, with the largest differences appearing for $\Omega_m$, $\omega_0$, and $\omega_a$. Within each catalog construction, however, the relative behavior of the cosmological approximations is similar: the Doublet reconstruction remains close to the exact result, the Singlet approximation shows larger departures in some redshift intervals, and the Fisher approximation generally produces broader credible intervals. Thus, although the GW-level approximation determines the information supplied to the cosmological analysis, the inclusion of second-order DALI terms at the cosmological level consistently reduces the discrepancy with respect to the exact likelihood.

The evolution of $\Delta H_0/H_0$ in Figure 11 isolates the effect of the different approximations on the Hubble constant constraint. To compare the complete cosmological parameter set at the final cumulative redshift, Table 5–Table 7 summarize the constraints obtained at $z_{\max}=1.0$. Focusing on $H_0$, the exact, Doublet, and Singlet approximations yield comparable statistical uncertainties, while the Fisher approximation consistently produces broader constraints across the three catalog configurations. The differences among the approximations are comparatively small in terms of $\Delta H_0/H_0$, indicating that the inferred precision on $H_0$ is relatively stable under these treatments. By contrast, the associated bias shows a stronger dependence on the adopted approximation and catalog construction. In particular, the Fisher approximation generally exhibits a smaller bias than the other methods, despite its larger uncertainty. These results indicate that the choice of approximation affects not only the width of the inferred $H_0$ constraints but also the displacement of the recovered value from the fiducial cosmology.

Method	$\Delta H_0/H_0$	$bias(H_0)$	$\Delta\Omega_M$	$bias(\Omega_M)$	$\Delta\omega_0$	$bias(\omega_0)$	$\Delta\omega_a$	$bias(\omega_a)$
Exact	3.40%	170.17%	0.092	107.92%	0.634	7.73%	4.290	80.92%
Doublet	3.50%	164.52%	0.097	76.28%	0.620	8.96%	4.468	55.51%
Singlet	3.36%	161.86%	0.174	39.26%	0.710	2.08%	5.655	39.04%
Fisher	5.90%	91.20%	0.760	16.48%	1.426	11.06%	19.781	16.66%

Method	$\Delta H_0/H_0$	$bias(H_0)$	$\Delta\Omega_m$	$bias(\Omega_m)$	$\Delta\omega_0$	$bias(\omega_0)$	$\Delta\omega_a$	$bias(\omega_a)$
Exact	2.54%	77.8%	0.08	123%	0.40	169%	3.07	172%
Doublet	2.69%	78.2%	0.07	152%	0.46	139%	3.29	159%
Singlet	3.31%	59.9%	0.09	84.3%	0.66	87.7%	5.66	53.3%
Fisher	5.42%	30.5%	0.19	64.1%	1.28	54.9%	12.51	58.4%

Jessen-Shannon Divergences Comparison

Across the three catalog constructions, the comparison shows a consistent ordering of the central JSD values, with the Doublet approximation yielding lower divergences than the Singlet and Fisher approximations for most cosmological parameters and redshift bins (Figure 12–Figure 14). The Singlet results generally occupy an intermediate range, while the Fisher approximation produces the largest divergences. This ordering is preserved despite differences in the redshift dependence among the cosmological parameters, indicating that the relative agreement of the three approximations with the exact posterior is primarily determined by the approximation scheme rather than by the catalog construction.

Method	$\Delta H_0/H_0$	$bias(H_0)$	$\Delta\Omega_M$	$bias(\Omega_M)$	$\Delta\omega_0$	$bias(\omega_0)$	$\Delta\omega_a$	$bias(\omega_a)$
Exact	4.82%	20.7%	0.14	130.7%	0.98	57.9%	8.91	66.4%
Doublet	5.16%	20.6%	0.17	88.0%	0.96	59.3%	10.72	14.2%
Singlet	4.60%	18.7%	0.15	104.1%	0.99	44.9%	13.23	5.0%
Fisher	19.23%	2.2%	0.67	33.0%	6.63	12.2%	101.74	3.5%

Method	$\Delta H_0/H_0$	$bias(H_0)$	$\Delta\Omega_M$	$bias(\Omega_M)$	$\Delta\omega_0$	$bias(\omega_0)$	$\Delta\omega_a$	$bias(\omega_a)$
Exact	2.54%	77.8%	0.08	123%	0.40	169%	3.07	172%
Doublet	2.69%	78.2%	0.07	152%	0.46	139%	3.29	159%
Singlet	3.31%	59.9%	0.09	84.3%	0.66	87.7%	5.66	53.3%
Fisher	5.42%	30.5%	0.19	64.1%	1.28	54.9%	12.51	58.4%

Evaluation Time

In addition to the accuracy of the different likelihood approximations, their computational cost is relevant for cosmological inference, where the likelihood must be evaluated repeatedly during posterior sampling. Figure 15 shows the distribution of the cosmological likelihood evaluation times for the exact, Doublet, and Singlet prescriptions, considering 80 catalogs containing 300 sources up to $z_{\max} = 1.0$.

The mean evaluation time obtained for the exact likelihood is $\bar{t}_{\mathrm{exact}} = 2.7 \times10^{-1}\ \mathrm{ms}$, compared with $\bar{t}_{\mathrm{Doublet}} = 2.3 \times10^{-2}\ \mathrm{ms}$ and $\bar{t}_{\mathrm{Singlet}} = 8.0 \times10^{-3}\ \mathrm{ms}$. The Doublet approximation is therefore approximately 12 times faster than the exact likelihood, while the Singlet approximation reduces the evaluation time by a factor of approximately 34. The distributions in Figure 15 also show a clear separation among the three prescriptions, indicating that the reduction in computational cost is systematic across the analyzed catalogs.

Although Singlet provides the lowest computational cost, the previous results show that the inclusion of the higher-order Doublet terms substantially improves the description of the non-Gaussian likelihood structure. Doublet therefore provides a favorable compromise in the present analysis: it retains part of the higher-order information beyond the Fisher approximation while remaining more than one order of magnitude faster to evaluate than the exact likelihood.

Conclusions and perspectives

In this paper, we have conducted a systematic analysis of the impact of different-quality bright siren simulations, using Gaussian and beyond-Gaussian likelihood approximations to constrain cosmological parameters. We found that, by using the best Hubble diagram approximation (via the GW-Doublet catalogue), the uncertainties in the cosmological parameters are reduced by approximately a factor of two, leading to tighter constraints. However, we also find that the bias in the cosmological parameters can be larger for the Doublet reconstruction, by factors of approximately $\sim2.4$–$2.7$ relative to the Fisher reconstruction, depending on the parameter. This is explained by the fact that the use of traditional priors in the GW analysis (particularly on $d_L$) can shift the measurements of $d_L$ away from the fiducial injections, which consequently increases the bias in the cosmological measurements. It is naturally expected that the use of Gaussian priors (centered on the injections) on GW parameters in the Doublet approximation can significantly reduce the bias in the cosmological parameters, leading to results similar to those obtained using Fisher with Gaussian priors. However, in this work, our goal was to investigate the performance of the DALI approximation as a computationally efficient representation of the traditional GW simulations based on exact likelihoods and traditional priors on the GW parameters.

Given the several GW approaches for Hubble diagrams, we also investigated the ability of the cosmological DALI formalism to recover results similar to those obtained using the exact cosmological likelihoods by comparing parameter uncertainties and Jensen–Shannon divergences. As expected, based on the Jensen–Shannon divergence, the Doublet approximation provides the closest approximation to the reference posterior among the methods considered. Within the Fisher-based approaches, Fisher sampling (Singlet-DALI) generally provides a closer approximation than the direct Fisher reconstruction. For the GW-Doublet catalogue, the cosmological uncertainties obtained with the Doublet likelihood differ from the exact results by only $\sim10\%$ on average across the four cosmological parameters, compared with $\sim195\%$ when using the traditional Fisher approximation. Therefore, the cosmological Doublet approximation is, on average, approximately 19 times closer to the exact likelihood results than the Fisher approximation.

Finally, we also investigated the computational cost of the different approximations at both the cosmological and GW levels. For the cosmological analysis, the Doublet likelihood evaluation is approximately 12 times faster than the exact likelihood, with mean evaluation times of $2.3 \times10^{-2}$ ms and $2.7 \times10^{-1}$ ms, respectively. To estimate the total computational cost of the combined GW and cosmological analyses, we further assume, as an approximate scaling, that the number of likelihood evaluations is proportional to the dimensionality of the parameter space. Since the GW analysis here involves six free parameters, compared with four for the cosmological analysis, we therefore take $N_{\mathrm{eval}}^{\mathrm{cosmo}}/N_{\mathrm{eval}}^{\mathrm{GW}} \simeq4/6$. Combining this assumption with the GW computational costs reported in previous work [50], where the Doublet approximation reduces the total GW MCMC cost by a factor of $\sim55$ relative to

the exact likelihood, we estimate that the total computational cost of the combined GW and cosmological analysis is reduced by approximately the same factor of $\sim55$. This is because the GW analysis dominates the total computational cost, while the cosmological Doublet provides an additional $\sim12\times$ reduction in the cost per likelihood evaluation. These estimates should be regarded as approximate, since the actual number of likelihood evaluations depends on the sampling algorithm and convergence criteria.

Overall, this work demonstrates the potential of DALI-based approximations to provide accurate and computationally efficient tools for forecasts in bright siren cosmology, while also highlighting the importance of carefully accounting for biases introduced by the underlying GW analysis. These results also motivate the exploration of similar approximation and optimization strategies in the context of dark sirens, where the cosmological likelihoods involve additional information from the galaxy distribution and can be substantially more computationally demanding. In this context, the use of efficient likelihood approximations and other computational optimization strategies may become particularly relevant for enabling cosmological inference with the large number of GW events expected from future detectors. Future work could extend this analysis to state-of-the-art waveform models for binary neutron stars that incorporate tidal effects and neutron-star tidal deformabilities, such as TEOBResumS [54] and waveform models incorporating the NRTidalv3 prescription [55], and investigate the impact of updated detector sensitivity curves on the resulting parameter-estimation and cosmological forecasts.

Acknowledgments

MHP is supported by the Brazilian agency Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). JMSdS is supported by the Brazilian research agency Fundação de Amparo à Pesquisa e Inovação do Estado do Espírito Santo (FAPES). He also thanks Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), process number E-26/200.236/2024 and E-26/200.237/2024. ODA thanks the Brazilian Ministry of Science, Technology and Innovation and the Brazilian Space Agency (AEB), which supported the present work under PO 20VB.0009. He also acknowledges CNPq for the Bolsa de Produtividade em Pesquisa (PQ), process number 310087/2021-0. We acknowledge the use of the computational resources of the joint CHE / Milliways cluster, supported by a FAPERJ grant E-26/210.130/2023. The authors also thank MCTIC/FINEP (CT-INFRA grant 011052700) and the Embrace Space Weather Program for the computing facilities at INPE.