# Dynamical dark energy models in the light of Gravitational-Wave Transient
Catalogues

Celia Escamilla-Rivera  0000-0002-8929-250X<EMAIL_ADDRESS>Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México,
Circuito Exterior C.U., A.P. 70-543, México D.F. 04510, México  Antonio Nájera
<EMAIL_ADDRESS>Instituto de Ciencias Nucleares,
Universidad Nacional Autónoma de México, Circuito Exterior C.U., A.P. 70-543,
México D.F. 04510, México

###### Abstract

The study of current gravitational waves catalogues provide an interesting
model independent way to understand further the nature of dark energy. Taking
advantage of them, in this work we present an update of the constraints
related to dynamical dark energy parameterisations using recent Gravitational-
Wave Transient catalogues (GWTC1 and GWTC-2). Also, we present a new treatment
for GW to establish the relation between the standard luminosity distance and
the siren distance. According to our Bayesian results developed with our join
SNeIa+CC+GW database, the $\Lambda$CDM model shows a preference against all
the dark energy parameterisations considered here. Moreover, with the current
GW transient database the GR standard luminosity and siren distances ratio
shows a strong preference against the modified gravity $\delta$-models
considered here.

###### pacs:

98.80.-k, 95.36.+x, 04.30.-w, 98.80.Es

##  I Introduction

One of the current challenges of the precision cosmology is to understand the
nature behind the late-cosmic acceleration. The study of a wide variety of
observations, e.g. Cosmic Microwave Background Radiation (CMBR), Supernovas
Type Ia (SNeIa), Baryonic Acoustic Oscillations (BAO) and the current
Gravitational Waves (GWs) have been useful to constraint the cosmological
parameters that define a specific cosmological model. Among these
observations, some of them are consistent with the standard $\Lambda$CDM,
which has its own theoretical and naturalness problems, so it is worthwhile to
consider alternative options. A convenient road is to keep General Relativity
(GR) as the gravitational theory and modify the stress-energy content of the
universe –the so-called dark energy– i.e. we modify the right-hand side of the
Einstein equations. A second road is to design an alternative theory of
gravity, whose additional degrees of freedom can handle the cosmic
acceleration –so-called the modified gravity method– via modifying the left-
hand side of the Einstein equations, which leads to new physics on small
cosmological scales. On the other hand, one of the most straightforward ways
to explore a specific cosmic accelerated expansion is at a phenomenological
level. In such a case, both roads can be quantified by the evolution of an
equation of state (EoS) parameterisation (geometric in the case of modified
gravity) of dark energy. Given the information on the values of the Hubble
parameter $H_{0}$, and the current matter density fraction, $\Omega_{m}$, an
arbitrary expression for the expansion evolution can be reproduced by assuming
a flat Friedmann-Robertson-Walker (FRW) space-time with a dark energy
component that has an EoS $w(z)$. Scenarios were a possible $w(z)\neq 1$ have
been discussed extensively and this denotes a strong evidence of new
gravitational physics and a possible new source of the cosmic acceleration.
Fitting a constant $w$ using observational data results in good agreement with
a EoS value of $-1$, but such analysis would have missed subtle variations in
the EoS, as for example if the average happened to be near to $-1$. Likewise,
using a specific parametric form of $w(z)$ is liable to biasing the outcome.
In this line of thought, some proposals have been currently studied, e.g.
bidimensional dark energy models (Escamilla-Rivera:2016qwv, ), EoS evolving
with a scalar field dark energy (Escamilla-Rivera:2016aca, ), generic EoS
derived from $f(R)$ (Jaime:2018ftn, ; Escamilla-Rivera:2020giy, ;
Cruz:2020cje, ), double exponential potentials for late-cosmic acceleration
(Das:2018xvj, ), interacting dark energy models like $w(q)$ (Elizalde:2018ahd,
), Padé approximants for dark energy (Rezaei:2017yyj, ) and dark energy
inverse cosmography (Escamilla-Rivera:2019aol, ; Munoz:2020gok, ).

Future proposals will be need to be addressed in the next decade by taking
into account the discordance present between the different cosmological probes
(DiValentino:2020vhf, ). Recently, constraints derived from considering the
interplay between GWs and other light gravitational degrees of freedom,
related to dark energy, have been used to strongly restrict the latter. Also,
with the help of the electromagnetic transient characteristic from the first
multi-messenger observations of a binary neutron star merger (GBM:2017lvd, ),
a possibility to constrain modified gravities emerged. In this matter, several
classes of scalar-tensor theory were found incompatible with the GW data or
else irrelevant to gravity at sufficiently large distances (Ezquiaga:2018btd,
). Some loopholes have been discussed by scanning a wider range of modified
gravity theories such as Beyond Horndeski and DHOST (Bordin:2020fww, ).

The goal of this work is to present an update in the study of dark energy
parameterisations by focusing on binary neutron stars with a counterpart
detected through a gamma-ray burst (GRB), using Gravitational-Wave Transient
Catalogues of Compact Binary Mergers observed by LIGO and Virgo during the
first and second observing runs (LIGOScientific:2018mvr, ) in order to shed
some light on their cosmological constraints.

This paper is organized as follows: In Sec.II we will explain how distances in
cosmology are written in standard way. This small introduction will encourage
the search for the theoretical siren distance definition. In Sec.II.5 we will
describe general properties to detect GWs and in Sec.III their standard
sirens. In Sec.IV we will present the theory behind the dark energy
parameterisations. Some alternative bidimensional parametrisations are
considered for the EoS, which deviate from the standard $\Lambda$CDM
$w(z)=-1$. Among them, we will study the Chevallier-Polarski-Linder (CPL)
(Chevallier:2000qy, ) and the Low Correlation (LC) (wang2008figure, )
parametrisations. In Sec.V we will the observational compilations considered
in this work. In Sec.VI, we will perform the statistical analysis over the
dark energy models discussed using the observational compilations including
the GW catalogues. Finally, in Sec.VII we will discuss our main results.

##  II Cosmological distances

###  II.1 Comoving distance

Comoving distance $D_{C}$ does not change with the expansion Hubble flow,
therefore travelling in a small time $dt$ as $dD_{C}=dD/a=cdt/a$
$\left(\text{since}\;c=dD/dt\right)$ we can express

$D_{C}(t)=c\int_{t}^{t_{0}}\frac{dt^{\prime}}{a(t^{\prime})},$ (1)

where $t$ is the time of light emission and $t_{0}$ is the time of detection.
This quantity can be also written as

$D_{C}(z)=c\int_{0}^{z}\frac{dz^{\prime}}{H(z^{\prime})},$ (2)

with $H(z)$ as the Hubble parameter $\left(H(t)=(1/a)da/dt\right)$. Defining
$E(z)\equiv H(z)/H_{0}$, (2) can be written as

$D_{C}(z)=D_{H}\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}.$ (3)

with $D_{H}$ the Hubble distance ($D_{H}=c/H_{0}$). As an example, the
dimensionless comoving distance $(D_{c}/D_{H})$ given for 4 different models
is described in Figure 1.

Figure 1: Examples of the dimensionless comoving distance. We present four
models: the $\Lambda$CDM model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.3,0.0,0.7)$, a positively curved
model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.6,-0.6,1.0)$, a
negatively curved model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.1,0.5,0.4)$ and the Einstein-de
Sitter model with $(\Omega_{m}\Omega_{K},\Omega_{\Lambda})=(1.0,0.0,0.0)$.

###  II.2 Transverse comoving distance

According to the $\Omega_{K}$ value, the transverse comoving distance
($D_{M}$) is defined as (distances, ):

$D_{M}(z)=\begin{cases}\dfrac{D_{H}}{(\Omega_{K})^{1/2}}\sinh\left[\dfrac{(\Omega_{K})^{1/2}D_{C}(z)}{D_{H}}\right],~{}\text{for}~{}\Omega_{K}>0,\\\
\\\ D_{C}(z),\quad\text{for}\quad\Omega_{K}=0,\\\ \\\
\dfrac{D_{H}}{(-\Omega_{K})^{1/2}}\sin\left[\dfrac{(-\Omega_{K})^{1/2}D_{C}(z)}{D_{H}}\right],~{}\text{for}~{}\Omega_{K}<0.\end{cases}$
(4)

The transverse comoving distance is defined as the ratio between the current
transverse velocity with its proper motion (weinberg1972gravitation, ). Notice
that if $\Omega_{K}=0$, we obtain an expression equivalent to the comoving
distance. In Figure 2 we discuss 4 possible cases for (4).

Figure 2: Examples of the dimensionless transverse comoving distance. We
present four models: the $\Lambda CDM$ model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.3,0.0,0.7)$, a positively curved
model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.6,-0.6,1.0)$, a
negatively curved model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.1,0.5,0.4)$ and the Einstein-de
Sitter model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(1.0,0.0,0.0)$.

###  II.3 Angular diameter distance

The angular diameter distance is the ratio between the physical size of an
object with its subtended angle, which in terms of the transverse comoving
distance is denoted by (dodelson2020modern, )

$D_{A}=\frac{D_{M}(z)}{1+z}.$ (5)

At small $z$ it is almost similar to the transverse comoving distance. The
dimensionless angular diameter distance $(D_{A}/D_{H})$ for our four
considered models is shown in Figure 3.

Figure 3: Examples of the dimensionless angular diameter distance. We present
four models: the $\Lambda$CDM model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.3,0.0,0.7)$, a positively curved
model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.6,-0.6,1.0)$, a
negatively curved model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.1,0.5,0.4)$ and the Einstein-de
Sitter model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(1.0,0.0,0.0)$.

###  II.4 Luminosity distance

Another way to compute the distance to a source is to measure its flux $F$ for
an object of known luminosity $L$. The luminosity distance $D_{L}$ is

$D_{L}=\left(\frac{L}{4\pi F}\right)^{1/2}.$ (6)

The luminosity distance can be written in terms of the angular diameter
distance and the transverse comoving distance as (dodelson2020modern, )

$D_{L}(z)=(1+z)D_{M}=(1+z)^{2}D_{A}.$ (7)

The dimensionless luminosity distance $(D_{L}/D_{H})$ for the four considered
models is shown in Figure 4.

Figure 4: Example of the dimensionless luminosity distance. We present four
models: the $\Lambda$CDM model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.3,0.0,0.7)$, a positively curved
model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.6,-0.6,1.0)$, a
negatively curved model with
$(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(0.1,0.5,0.4)$ and the Einstein-de
Sitter model with $(\Omega_{m},\Omega_{K},\Omega_{\Lambda})=(1.0,0.0,0.0)$.

###  II.5 Standard sirens distance

In GR, GWs have two transverse polarizations (schutz2009first, ). Each
polarization are denoted by $h_{+}$ and $h_{\times}$.

Let us design a ring of free particles where a GW arrive. By considering two
well-chosen test masses (see Figures 5 and 6) we can calculate their distance
to the center to see if a GW is causing a perturbation. In the first half
cycle, the distance to the center of one test mass will increase, while the
other will decrease. In the second half cycle, it is on reverse. By building
an interferometer a passing space-time perturbation would produce a time-
varying distance that can be analysed to detect the GW. Unfortunately, the
changes in distances by astrophysical situations are minimal because of the
tiny amplitudes (dhurandhar2017cosmic, ). This is the reason of why a very
sensitive interferometer is required in the experiment.

Figure 5: (a) A ring of particles initially at rest. (b) The ring of particles
distortioned by a $h_{+}$ at $\pi/4$ phase. (c) The distortion at phase
$3\pi/2$. Two test masses are being considered.

If a GW with $h_{+}$ that is traveling in the z-direction hits a circular ring
initially at rest (the phase of the GW when it hits the ring is assumed 0) in
the $x-y$ plane the particles will be moved in the way shown in Figure 5.
Since the GW is periodic, (b) and (c) will be repeated every half cycle. The
case for a $h_{\times}$ wave is shown in Figure 6. An arbitrary GW will be a
linear combination of these two polarization states.

Figure 6: (a) A ring of particles initially at rest. (b) The ring of particles
distortioned by a $h_{\times}$ at $\pi/4$ phase. (c) The distortion at phase
$3\pi/2$. Two test masses are being considered.

The tensor $h_{ij}$ of a GW can be expressed (this derivation is performed in
Appendix A) as

$h_{ij}=\frac{2G}{c^{4}D_{S}}\left[\frac{d^{2}I_{ij}\left(t-\dfrac{R}{c}\right)}{dt^{2}}\right]^{TT},$
(8)

where $I_{ij}$ denotes the quadrupole moment tensor and $TT$ is the transverse
and traceless part of $I_{ij}$. By integrating the instantanous energy flux
gravitational radiation over all directions the luminosity of a GW can be
computed as (damour2016gravitational, )

$\mathcal{L}=\frac{G}{5c^{2}}\left[\frac{d^{3}I_{ij}^{TT}}{dt^{3}}\frac{d^{3}I^{ij\;TT}}{dt^{3}}\right].$
(9)

Hence the time average over a period can be written as

$\mathcal{L}=\frac{G}{5c^{2}}\left<\frac{d^{3}I_{ij}^{TT}}{dt^{3}}\frac{d^{3}I^{ij\;TT}}{dt^{3}}\right>.$
(10)

The transverse and traceless part was taken into account, which means that it
is being considered the non-spherical motion of the GW source, hence only the
kinetic energy in the non-spherical motion contributes to the emission of
energy of a GW.

We can derive the amplitude of the system considering an inspiralling binary
system with $M$ as the total mass. Also, it is convenient to consider their
reduced mass $\mu=m_{1}m_{2}/M$, since the kinetic energy of the system is
$K=\mu v^{2}/2$. Using Kepler’s third law we can compute the amplitude of a GW

$\frac{T^{2}}{r^{3}}=\frac{4\pi^{2}}{GM},$ (11)

where $T$ and $r$ denotes the period of the orbit and the distance between the
masses, respectively. We can rewrite the latter equation in terms of the
velocity of rotation if $v=2\pi r/T$ and $v^{2}=GM/r$. By inserting the
kinetic energy of the system into (8) the amplitude of a GW is
(dhurandhar2017cosmic, )

$h\approx\frac{G^{2}\mu M}{D_{S}rc^{4}},$ (12)

where $D_{S}$ is the siren distance to the source and $r$ the distance between
the binaries. As an example, we can consider a couple of $10M_{\odot}$ mass
binaries, so that the total mass is $20M_{\odot}$ at a siren distance of $200$
Mpc and separated by $r=GM_{\odot}/c^{2}$, therefore the amplitude is
$h\approx 2.39\,\times\,10^{-20}$.

##  III Dynamics

Using (10), the luminosity of a GW created by an inspiralling binary is
(damour2016gravitational, )

$\mathcal{L}=\frac{32\mu^{2}G^{4}M^{3}}{5r^{5}c^{5}}.$ (13)

This is the observational luminosity, from which the observational luminosity
siren distance can be computed by deriving the dynamics of the system. Since
binaries emit energy in the form of GWs, their kinetic energy must be reduced
and therefore the orbit radius decreases with time, and eventually the
binaries collide with each other. To compute this rate of change of the orbit
radius we can consider the following ansatz: if the entire bond energy of the
system (which is $E=\mu v^{2}/2$) transforms into luminosity energy, then

$\mathcal{L}=-\frac{dE}{dt}.$ (14)

Using Kepler’s third law to write the bond energy in function of the radius
$E=\mu GM/2r$, where $M$ is the total mass of the system, therefore

$\mathcal{L}=-\frac{\mu GM}{r^{2}}\frac{dr}{dt}.$ (15)

Solving for $dr/dt$ we obtain

$\frac{dr}{dt}=-\frac{64G^{3}\mu M^{2}}{5r^{3}c^{5}},$ (16)

which can be rewritten in terms of the angular frequency of rotation and
$T=2\pi/\omega$, the distance between the binaries is given by

$r=\left(\frac{GM}{\omega^{2}}\right)^{1/3},$ (17)

and (16) in terms of the angular frequency is

$\frac{d\omega}{dt}=\frac{96\mu}{5c^{5}}\left(\omega^{11}G^{5}M^{2}\right)^{1/3}.$
(18)

The orbital and GW frequency of the binaries are related by
$f_{\text{GW}}=f_{\text{orbital}}/2=\omega/\pi$, then

$\frac{df_{\text{GW}}}{dt}=\frac{96\mu\pi^{8/3}}{5c^{5}}\left(f_{\text{GW}}^{11}G^{5}M^{2}\right)^{1/3}.$
(19)

Combining (12) and (19) we can obtain the siren distance for a GW

$D_{S}=\frac{5c}{96\pi^{2}hf_{\text{GW}}^{3}}\frac{df_{\text{GW}}}{dt}.$ (20)

This expression allows us to measure the siren luminosity distance, the
frequency and the strain as a function of time111From this we can determine
the quantity $df_{\text{GW}}/dt$ (abbott2016observation, ).. The amplitude $h$
and the strain are related by solving the geodesic deviation equation, then
$h\approx\Delta L/L$ (where $L$ is the distance between some test masses and
$\Delta L$ is the strain). Therefore, Eq.(20) can be computed with the
observables of a GW. The time it takes to the binaries to coalesce can also be
computed by integrating (19). The lifetime of the system (so called the chirp
time ($t_{c}$)) from an initial radius $r_{0}$ at an initial frequency
$f_{0\;\text{GW}}$ is given by

$t_{c}=\frac{5c^{5}}{256\mu}\left(\frac{1}{\pi^{8}f_{0\;\text{GW}}^{8}G^{5}M^{2}}\right)^{1/3}.$
(21)

In comparison to the other definitions of distances described, a theoretical
siren luminosity distance is required at this point. The experimental distance
was derived in (20) and experiments that detect GWs report them along with the
redshift of the event, so that the theoretical is the quantity to determine. A
possible form is proposed in (Mitra:2020vzq, ) where a relationship between
the standard luminosity and the siren luminosity distance is given by

$D_{S}(z)=D_{L}(z)\exp\left(-\int_{0}^{z}\frac{\delta(z^{\prime})}{1+z^{\prime}}dz^{\prime}\right),$
(22)

which comes from a modified gravity model. The subscript $S$ and $L$ denotes
the siren luminosity distance and the standard candle luminosity distance,
respectively. To analyse data coming from GWs, the form of $\delta(z)$ needs
to be found.

In this work we consider three forms of delta functions from
(belgacem2019testing, ). These deltas are defined as:

$\delta_{1}(z)=\frac{n(1-\Xi_{0})}{1-\Xi_{0}+\Xi_{0}(1+z)^{n}},$ (23)
$\delta_{2}(z)=\frac{n(1-\Xi_{0})(1+z)z^{n-1}}{1+\Xi_{0}(\exp(z^{n})-1)},$
(24)
$\delta_{3}(z)=\frac{n(1-\Xi_{0})}{1-\Xi_{0}+\Xi_{0}(1+z)^{n}}+\dfrac{n(1-\Xi_{0})}{(1+z)^{n}}-\dfrac{2n(1-\Xi_{0})}{(1+z)^{2n}},$
(25)

where the quantities $n$ and $\Xi_{0}$ are both positive and given by
(Mitra:2020vzq, )

$\begin{split}\Xi_{0}=\lim_{z\to\infty}\frac{M^{*}(0)}{M^{*}(z)},\\\
n\approx\frac{\alpha_{M0}}{2(\Xi_{0}-1)}.\end{split}$ (26)

The GR case can be recover when $\Xi_{0}=1$ (belgacem2019testing, ) and its
delta is

$\delta_{\text{GR}}(z)=0.$ (27)

By inserting (23), (24), and (25) into (22) and integrating, the siren
distances are

$D_{S_{1}}(z)=D_{L}(z)\left(\Xi_{0}+\dfrac{1-\Xi_{0}}{1+z^{n}}\right),$ (28)
$D_{S_{2}}(z)=D_{L}(z)\left[\Xi_{0}+(1-\Xi_{0})\exp(-z^{n})\right],$ (29)
$\displaystyle
D_{S_{3}}(z)=D_{L}(z)\exp\left(-\dfrac{(1-\Xi_{0})(1-(1+z)^{n})}{(1+z)^{2n}}\right)\times\left(\Xi_{0}+\dfrac{1-\Xi_{0}}{1+z^{n}}\right),$
(30)

Finally, the GR case can be achieved by the ratio

$D_{S_{\text{GR}}}(z)=D_{L}(z).$ (31)

It is useful to analyse the limits of (28), (29), and (30): when $z\to 0$,
$D_{S}(z)=D_{L}(z)$, and when $z\to\infty$, $D_{S}(z)\to D_{L}(z)\,\Xi_{0}$.
Therefore, at $z=0$, both distances are equal and at early times their ratio
goes to a positive constant. Also, if $\Xi_{0}=1$, $D_{S}(z)=D_{L}(z)$, the
ratio with these distances holds for every $z$. In this work we also develop
these parameterisations over $\delta$ to study the consequences in the dark
energy dynamics.

##  IV Dynamical dark energy parameterisations

Since the concordance cosmological $\Lambda$CDM model has several theoretical
problems (Weinberg:2000yb, ), new proposals have been presented along these
years. One of the most popular is related to the addition of an exotic fluid
denoted as dark energy. Furthermore, modifications of the gravity itself
(Clifton:2011jh, ) also have been considered.

Several dark energy EoS $w$ have been built to describe the late cosmic
acceleration. To see which one of them has more constraining power we test
them using observational data.

The fact that the universe is experimenting an accelerated expansion direct us
to consider a negative pressure component at late times. Since the energy is
always positive, the dark energy EoS $w=w(z)$ requires to be
$(w(z)=P/\rho<-1/3)$. For a spatially flat universe

$E(z)=\left(\frac{H(z)}{H_{0}}\right)=\left[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}\,f(z)\right]^{1/2},$
(32)

where the function $f(z)$ is defined as

$f(z)=\exp\left(3\int_{0}^{z}\frac{1+w(z^{\prime})}{1+z^{\prime}}dz^{\prime}\right),$
(33)

for a given $w(z)$ the dynamics of the expansion is determined. The
$\Lambda$CDM model is a special case when $w(z)=-1$, $f(z)=1$. From here, we
can consider the following six bidimensional222We denoted bidimensional as the
dependence solely of two free theoretical parameters $w_{i}$ dark energy
parameterisations 333A previous analysis of these models have been presented
in (Escamilla-Rivera:2016qwv, ) with different analytical forms of $w(z)$.

###  IV.1 Lambda Cold Dark Matter-redshift $(\mathbf{\Lambda CDM})$

This cosmological model is represented by $w(z)=-1$, so the Friedmann equation
becomes

$E(z)^{2}=\Omega_{m}(1+z)^{3}+(1-\Omega_{m}).$ (34)

This is a simple model that describes the late-time cosmic acceleration.

###  IV.2 Chevallier-Polarski-Linder ($\mathbf{CPL}$)

Another parameterization considers an EoS (Chevallier:2000qy, ; Linder:2007wa,
) of the form

$w(z)=w_{0}+\left(\frac{z}{1+z}\right)w_{1},$ (35)

which converges for high redshift and also reduces to the $\Lambda CDM$ model
for $w_{1}=0$ and $w_{0}=-1$. In this case the Friedmann equation is

$\displaystyle E(z)^{2}$ $\displaystyle=$
$\displaystyle\Omega_{m}(1+z)^{3}+(1-\Omega_{m})\exp\left(-\frac{3w_{1}z}{1+z}\right)$
(36) $\displaystyle\times(1+z)^{3(1+w_{0}+w_{1})}.$

###  IV.3 Barboza-Alcaniz ($\mathbf{BA}$)

This model proposes an EoS (barboza2008parametric, ) of the form

$w(z)=w_{0}+\frac{z(1+z)}{1+z^{2}}w_{1}.$ (37)

Notice that this expression does not diverge for $z\to-1$. In this case the
Friedmann equation is

$\displaystyle E(z)^{2}$ $\displaystyle=$
$\displaystyle\Omega_{m}(1+z)^{3}+(1-\Omega_{m})(1+z)^{3(1+w_{0})}$ (38)
$\displaystyle\times(1+z^{2})^{3w_{1}/2}.$

###  IV.4 Jassal-Bagla-Padmanabhan ($\mathbf{JBP}$)

This parameterization has an EoS given by (jassal2005wmap, )

$w(z)=w_{0}+\frac{z}{(1+z)^{2}}w_{1},$ (39)

which converges for high redshifts and contains a dark energy contribution at
this scale. The Friedmann equation for this parameterization is

$\displaystyle E(z)^{2}$ $\displaystyle=$
$\displaystyle\Omega_{M}(1+z)^{3}+(1-\Omega_{M})$ (40)
$\displaystyle\times\exp\left(\frac{3w_{1}z^{2}}{2(1+z)^{2}}\right)(1+z)^{3(1+w_{0})}.$

###  IV.5 Linear-redshift ($\mathbf{LR}$)

The EoS for this model is given by $w(z)=w_{0}-w_{1}z$ (PhysRevD.64.123527, )
and it can be reduced to $\Lambda$CDM when $w_{0}=-1$ and $w_{1}=0$. With this
EoS the $f(z)$ can be computed and the Friedmann equation becomes

$E(z)^{2}=\Omega_{m}(1+z)^{3}+(1-\Omega_{m})\exp(-3w_{1}z)(1+z)^{3(1+w_{0}+w_{1})}.$
(41)

For high redshift this EoS diverges. This gives strong constraints on the
value of $w_{1}$.

###  IV.6 Low Correlation ($\mathbf{LC}$)

This is a model with an EoS (wang2008figure, ) given by

$w(z)=\frac{(-z+z_{c})w_{0}+z(1+z_{c})w_{c}}{(1+z)z_{c}},$ (42)

where $w_{0}=w(0)$ and $w_{c}=w(z_{c})$ and $z_{c}$ is a given redshift value
so that $(w_{0},w_{z})$ are uncorrelated. These values depend on the data used
to be constrained. This parameterization reduces to $\Lambda$CDM for
$w_{0}=-1$ and $w_{c}=-1$ .The Friedmann equation is

$\displaystyle E(z)^{2}=\Omega_{m}(1+z)^{3}+(1-\Omega_{m})$
$\displaystyle\times\exp\left(\frac{9(w_{0}-w_{c})z}{1+z}\right)(1+z)^{3(1-2w_{0}+3w_{c})}.$
(43)

###  IV.7 Wetterich-redshift ($\mathbf{WP}$)

This parameterization is given by (wetterich2004phenomenological, )

$w(z)=\frac{w_{0}}{(1+w_{1}\ln(1+z))^{2}},$ (44)

where the Friedmann equation is

$\displaystyle E(z)^{2}$ $\displaystyle=$
$\displaystyle\Omega_{m}(1+z)^{3}+(1-\Omega_{m})$ (45)
$\displaystyle\times(1+z)^{3\left(1+\dfrac{w_{0}}{1+w_{1}\ln(1+z)}\right)}.$

##  V Observational compilations

It is well known that the dark energy domination era began at a late period in
the evolution of the universe. Then the best data to study this dark component
is at low redshifts. In this section we explain the method to obtain the
cosmological parameters $\Omega_{m}$, $\Theta=\\{w_{0},w_{1}\\}$,
$h=H_{0}/(100\,\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1})$. We perform these
analyses using the $\Lambda$CDM model and the six bidimensional models
described. We start by using the SNeIa and Cosmic Chronometers compilations,
once this is done, we will implement the GW catalogues data to improve the
constraints over the cosmological parameters.

###  V.1 Pantheon supernovae compilation

The Pantheon catalogue (scolnic2018complete, ) includes
$N_{\text{Pantheon}}=1048$ events. These events are treated in an interval
$z\in[0.010,2.26]$. To compute the statistical analysis of this sample we use

$\mu(z_{i},\mu_{0})=5\log_{10}\left(\dfrac{D_{L}(z)}{10\,\text{pc}}\right),$
(46)

where $D_{L}(z)$ is is given by (7). The best fit parameters $\Omega_{m}$,
$\Theta=\\{w_{0},w_{1}\\}$, $h$ can be calculated by maximizing the logarithm
of the likelihood function (hogg2010data, )

$\displaystyle\ln
p(\mu(z)|z,\sigma,h,\Omega_{m},\Theta)=-\frac{1}{2}\left(\chi^{2}+\sum_{n=1}^{N_{\text{Pantheon}}}\ln(2\pi\sigma_{\text{obs}}^{2})\right),$
(47)

where the $\chi^{2}$ is given by

$\chi^{2}=\sum_{n=1}^{N_{\text{Pantheon}}}\frac{(\mu(z_{n},\Omega_{m},\Theta)-\mu_{\text{obs}}(z_{n}))^{2}}{\sigma_{\text{obs}}^{2}},$
(48)

with $\sigma_{n\,\text{obs}}^{2}$ the measurement variances.

###  V.2 Cosmic Chronometers compilation

This database includes $N_{\text{CC}}=52$ events in the redshift interval
$z\in[0,2.36]$ (10.1093/mnras/sty260, ). In comparison to SNeIa data, the
logarithm of the likelihood function (47) is expressed as the difference
between the function $H(z)$ instead of $\mu(z)$. The $\chi^{2}$ function for
this case given by

$\chi^{2}=\sum_{i=1}^{N_{\text{Hubble}}}\frac{(H(z_{i},\Omega_{M},w_{0},w_{1})-H_{\text{obs}}(z_{i}))^{2}}{\sigma_{n\,\text{obs}}^{2}},$
(49)

where and $H_{\text{obs}}$ and $\sigma_{n\,\text{obs}}$ correspond to the
observations and their errors, respectively.

###  V.3 Gravitational-Wave Transient Catalogue GWTC-1

This catalogue includes $N_{\text{GWTC-1}}=11$ confident detection events in
redshift range of $z\in[0.01,0.49]$ from gravitational-wave searches for
coalescing compact binaries, all of them with masses $M>1M_{\odot}$. These
events observed by LIGO 444https://www.ligo.caltech.edu and
Virgo555http://www.virgo.infn.it come from the first and second observing runs
(abbott2019gwtc, ). In Table III from (abbott2019gwtc, ), the 11 events are
reported along with their observational siren distances $D_{S}$ in Mpc and
redshifts $z$ and other several observational variables, all of them with
their corresponding $90\%$ credible variances.

With the latter information, the natural logarithm of the likelihood function
for the GWTC-1 catalogues can be compute using the following expression

$\displaystyle\ln
p(D_{S\,\text{GWTC-1}}|z,\sigma,h,\Omega_{m},w_{0},w_{1},n,\Xi_{0})=-\frac{1}{2}\left(\chi^{2}+\displaystyle\sum_{j=1}^{N_{\text{GWTC-1}}}\ln(2\pi\sigma_{j\,\text{obs,
GWTC-1}}^{2})\right),$ (50)

where

$\chi_{\text{GWTC-1}}^{2}=\sum_{j=1}^{N_{\text{GWTC-1}}}\frac{(D_{S}(z_{j},h,\Omega_{M},w_{0},w_{1},n,\Xi_{0})-D_{S\,\text{obs,
GWTC-1}}(z_{j}))^{2}}{\sigma_{j\,\text{obs, GWTC-1}}^{2}}.$ (51)

The GWTC-1 catalogue includes variances in $z$, therefore the squared
variances are

$\sigma_{j\,\text{obs, GWTC-1}}^{2}=\sigma_{j\,D_{S}\text{obs,
GWTC-1}}^{2}+\sigma_{j\,z\text{obs, GWTC-1}}^{2}.$ (52)

Notice that the squared variances of the siren distance per redshift are
added. But since the upper and lower reported variances are not equal, the
biggest of each case were taken as the variance. Then the logarithm of the
likelihood function is maximized and we run Markov Chains in order to get the
best fit values of $h$, $\Omega_{M}$, $w_{0}$, $w_{1}$, $\Xi_{0}$, and $n$ and
their confidence regions.

###  V.4 Gravitational-Wave Transient Catalogue GWTC-2

This second GW catalogue includes $N_{\text{GWTC-2}}=39$ candidate GW events,
with less than $10\%$ of contamination fraction in the redshift range
$z\in[0.03,0.8]$ from the first half of the third observable run
(Abbott:2020niy, ). In Table VI from (Abbott:2020niy, ) these events measured
by LIGO-Virgo are presented with their observational siren distances $D_{S}$
in Gpc666Our units conversion must be performed to be consistent with the
GWTC-1 database, where the $D_{S}$ is reported in Mpc. and redshifts $z$,
along with other several observational variables, all of them with their
corresponding $90\%$ credible variances.

The logarithm of the likelihood function for the GWTC-2 catalogue is given by

$\displaystyle\ln
p(D_{S\,\text{GWTC-2}}|z,\sigma,h,\Omega_{M},w_{0},w_{1},n,\Xi_{0})=-\frac{1}{2}\left(\chi^{2}+\displaystyle\sum_{j=1}^{N_{\text{GWTC-2}}}\ln(2\pi\sigma_{j\,\text{obs,
GWTC-2}}^{2})\right),$ (53)

where

$\chi_{\text{GWTC-2}}^{2}=\sum_{j=1}^{N_{\text{GWTC-2}}}\frac{(D_{S}(z_{j},h,\Omega_{m},w_{0},w_{1},n,\Xi_{0})-D_{S\,\text{obs,
GWTC-2}}(z_{j}))^{2}}{\sigma_{j\,\text{obs, GWTC-2}}^{2}}.$ (54)

The GWTC-2 catalogue includes variances in $z$ described by

$\sigma_{j\,\text{obs, GWTC-2}}^{2}=\sigma_{j\,D_{S}\text{obs,
GWTC-2}}^{2}+\sigma_{j\,z\text{obs, GWTC-2}}^{2}$ (55)

Both squared variances of the siren distance and redshift are considered. And
as in the GWTC-1 database the upper and lower reported variances are not
equal, then the biggest of each case was taken as the variance. Then logarithm
of the likelihood function is maximized and we run again Markov Chains in
order to get the best fit values of $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$,
$\Xi_{0}$, and $n$ and their confidence regions. The analysis is done with the
sum of the GWTC-1 and GWTC-2 databases, we denote the join likelihood as ‘GW
database’. Due that this database reports events with less than $10\%$ of
contamination fraction, the best fit values and their confidence regions will
be biased a little. Moreover, this bias will be minimal due the density of
data when combining it with Pantheon and CC data, so that the GWTC-2 that
contains the bias represents just the $3.39\%$ of the total Pantheon+CC+GW
database.

##  VI Methodology

Once the best cosmological parameters were obtained by maximizing the
logarithm of the likelihood function, the Monte Carlo Markov chains to get the
$\sigma$ and $2\sigma$ confidence regions were obtained using a modified
version of emcee (foreman2013emcee, ). The confidence regions plots and best
fit values with their 1-$\sigma$ uncertainties were made by using the GetDist
software ((Lewis:2019xzd, ) and references therein).

###  VI.1 Statistical analysis

The best fit values and the 1-$\sigma$ uncertainties for the $h$,
$\Omega_{m}$, $w_{0}$, and $w_{1}$ parameters using Pantheon and CC databases
and their overlapping are reported in Table 2. The best fit values and the
1-$\sigma$ uncertainties for the $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$,
$\Xi_{0}$, $n$ parameters for the case of the GW database along with their
overlapping with Pantheon and CC are reported in Tables 3, 4 and 5, one for
each siren delta function $\delta_{1},\delta_{2}$ and $\delta_{3}$,
respectively.

We also include in the analysis the plots at 1-$\sigma$ and 2-$\sigma$
confidence regions for the seven models (including $\Lambda$CDM model with the
six bidimensional dark energy parameterisations). These analyses include for
each model:

  1. 1.

Analysis performed using Pantheon compilation, CC data, and their overlap
Pantheon+CC. See Figures 7-9-11-13-15-17-19.

  2. 2.

Analysis performed per siren delta function using Pantheon compilation, the
overlap Pantheon+CC and GW data. See Figures 8-10-12-14-16-18-20.

  3. 3.

Analysis performed per siren delta function using GW data, the overlap
Pantheon+GW, and the overlap Pantheon+CC+GW. See Figure 21.

###  VI.2 Information Criteria analysis

We performed a Markov-chain Monte Carlo (MCMC) analysis using a modified
version of the emcee code777https://emcee.readthedocs.io interfaced to the
publicly available sampling code
MontePython-v3888https://github.com/brinckmann/montepython_public
(Audren:2012wb, ; Brinckmann:2018cvx, ). We obtained plots using the GetDist
package999https://getdist.readthedocs.io/en/latest (Lewis:2019xzd, ). We
sampled over the parameters $\\{h_{0},\Omega_{m},\Theta,\Xi_{0},n\\}$ using
the Metropolis-Hastings method and with a Gelman-Rubin (Gelman:1992zz, )
convergence criteria around $R-1<0.03$. To quantify to what extent the
improvement in our fit to the database described warrants the increase in the
model complexity compared to the $\Lambda$CDM model, we compute the Bayes
factor –the ratio of the evidences for the extended model ${\cal M}_{\rm E}$
with respect to the baseline ${\cal M}_{\rm L}$– as (Jeffreys:1939xee, ):

$\displaystyle B_{EL}\equiv\frac{\int
d\bm{\theta}_{E}\,\pi(\bm{\theta}_{E}|{\cal M}_{E}){\cal
L}(\mathbf{x}|\bm{\theta}_{E},{\cal M}_{E})\,,}{\int
d\bm{\theta}_{L}\,\pi(\bm{\theta}_{L}|{\cal M}_{L}){\cal
L}(\mathbf{x}|\bm{\theta}_{L},{\cal M}_{L})\,,}\,,$ (56)

where $\pi(\bm{\theta}_{\rm E,\,L})$ is the prior for the parameters
${\theta}_{\rm E,\,L}$ and ${\cal L}(\mathbf{x}|\bm{\theta}_{\rm E,\,L})$ the
likelihood of the data given the model ${\cal M}_{\rm E,\,L}$. The extent to
what the extended model ${\cal M}_{\rm E}$ is preferred over the baseline
${\cal M}_{\rm L}$ can be qualitatively assessed using the Jeffreys scale
reported in Table 1. We compute the evidence directly from our MCMC using the
method introduced in Ref. (Heavens:2017afc, ) implemented in the MCEvidence
code101010https://github.com/yabebalFantaye/MCEvidence.

We compared the six dark energy parameterisations against the $\Lambda$CDM
model with the Pantheon+CC+GW dataset using the Akaike Information Criteria
(AIC) and the Bayesian Information Criteria (BIC) (Liddle:2004nh, ;
doi:10.1080/01621459.1995.10476572, ). The AIC is given by (Liddle:2004nh, )

$\text{AIC}=\chi_{\text{Pantheon+CC+GW min}}^{2}+2\,k,$ (57)

where $\chi_{\text{Pantheon+CC+GW min}}^{2}$ is the minimum $\chi^{2}$ given
by the sum of $\chi^{2}$s of the Pantheon+CC+GW data, and $k$ the number of
parameters of the model. Meanwhile, the BIC is given by (Liddle:2004nh, )

$\text{BIC}=\chi_{\text{Pantheon+CC+GW min}}^{2}+k\,\ln N,$ (58)

where $N$ is the number of elements of the dataset, for Pantheon+CC+GW it is
$N=1150$. When comparing models, the one with smaller AIC and BIC is preferred
and the strenght of the preference is given by the Jeffrey’s criterion (see
Table 1) (doi:10.1080/01621459.1995.10476572, ; belgacem2019testing, ) with
$\Delta\text{AIC}=|\text{AIC}_{\text{model 1}}-\text{AIC}_{\text{model 2}}|$
and $\Delta\text{BIC}=|\text{BIC}_{\text{model 1}}-\text{BIC}_{\text{model
2}}|$. We present the results of these analyses in Tables 6, 7 and 8.

Information Criteria | Preference for model ${\cal M}_{i}$
---|---
$0<\Delta\text{BIC},~{}\Delta\text{AIC}<2\qquad$ | Weak
$2\leq\Delta\text{BIC},~{}\Delta\text{AIC}<6$ | Definite
$6\leq\Delta\text{BIC},~{}\Delta\text{AIC}<10$ | Strong
$10\leq\Delta\text{BIC},~{}\Delta\text{AIC}<12\qquad$ | Very strong
Table 1: Revised Jeffreys scale used to interpret the values of $\ln B$
obtained when comparing two competing models through their Bayesian evidence.
A value of $\ln B>0$ indicates that the extended model is favoured with
respect to the $\Lambda$CDM baseline model.

To see whether one of the three siren delta functions (23), (24), (25) or the
ratio between the siren distance and luminosity distance ($D_{S}(z)=D_{L}(z)$)
derived from GR is better, we performed the AIC and BIC analyses by comparing
the three deltas against GR. For this analysis we considered the constrains
obtained from the six dark energy models along with the $\Lambda CDM$ model
using the full GW database. We present these results in Table 9.

As it can be seen in the plots, the $\Lambda$CDM model ($w_{0}=-1$) and the LC
model ($w_{1}=0$, $w_{c}=-1$) lies within the 2-$\sigma$ uncertainties for all
the cases and databases. Also, for all the cases and databases containing GW
data, the equivalence of the siren distance $D_{S}$ with the luminosity
distance $D_{L}$ ($D_{S}=D_{L}$) lies within the 1-$\sigma$ uncertainty.
Therefore, this ratio is consistent with the current GW data, and also it is
the best model.

According to the values reported in Table 3, current GW transient data is not
able to constrain the value of $\Omega_{m}$, $w_{1}$ and $n$. When an overlap
of GW with other database was done the values of $\Omega_{m}$, $w_{1}$ were
almost GW-non-dependent. Also, since the definitions for Pantheon and CC data
do not contain the $n$ parameter, a best fit value for this parameter cannot
be computed, because the form of the graphic of $n$ in all cases was non-
gaussian like. However, GW data was able to constrain the value of $\Xi_{0}$,
which is consistent with the ratio $D_{S}=D_{L}$, and gives a clue that both
distances are equivalent. Future GW databases will decide either this ratio
holds or not, and also might constrain the $n$ value.

As we can see in Tables 6, 7, and 8, since $\Lambda$CDM has the smallest AIC
and BIC in all cases (for the three deltas and the six dark energy models), it
is the best model with the current Pantheon+CC+GW data. How preferred is the
$\Lambda$CDM model over the six dark energy parameterisations was computed
with the Jeffrey’s criterion. For the AIC criterion the preference is very
weak, but with the BIC criterion $\Lambda$CDM has a very strong preference
against the other models. Hence, the best model by far with our dataset is
$\Lambda$CDM. But, among the dark energy parameterisations, BA was the
preferred model because it has the lower AIC and BIC values, for the three
deltas.

In the case of the deltas against GR, their values are reported in Table 9, GR
is the best model ($D_{L}(z)=D_{S}(z)$) since it has in all the models the
smallest AIC and BIC. The GW database is small ($N_{\text{GW}}=50$) and the
GTC-2 catalogues that forms part of this database include a contamination of
less than $10\%$, so that a bias of less than $10\%$ must be present in the
analysis. Therefore, the results should not be extremely accurate, but with
AIC, GR has a positive preference over $\delta_{1}$, $\delta_{2}$, and
$\delta_{3}$. Meanwhile, according to the BIC, GR has a strong preference over
the same deltas. Therefore the GR case is the best with the current GW
database. Furthermore, the luminosity distance ($D_{L}$) and the siren
distance ($D_{S}$) ratio is better than the expressions (28), (29), and (30)
with the current GW database. Among the three deltas, $\delta_{3}$ was the
best by a minimum difference because it has the lowest AIC and BIC, except for
the CPL case in which $\delta_{2}$ is the best parameterisation.

Furthermore, we computed the ratio between the siren distance and the
luminosity distance for the case $\delta_{3}$ using the AIC and BIC criteria
for the $\Lambda$CDM model, the six dark energy parameterisations, and the GR
case. This ratio is shown in Figure 22 using the delta functions for these
models.

##  VII Conclusions

Our first task in this work was to present an update analysis of the six dark
energy parameterisations using the new GW databases described. We consider
current SNeIa (Pantheon)+CC+GW data, the GW data came from the GWTC-1
(abbott2019gwtc, ) and GWTC-2 (Abbott:2020niy, ) catalogues. With this new
database, the $\Lambda$CDM had a preference against the six dark energy
considered models according to the BIC.

The second task was to study the siren distance of GWs. This was performed
using the siren delta functions that give three different siren distance -
luminosity distance relationships. They were compared against the GR ratio
case, which had a positive preference with the AIC and strong preference with
the BIC. As a result, the ratio between the siren distance and the luminosity
distance is the best model.

Although the standard models ($\Lambda$CDM and GR) were preferred against the
bidimensional dark energy parameterisations and the ratio between siren
distance-luminosity distance, GWs have opened the door to update the previous
analysis with supernovae and CC data with GW data. Future GW data might give
more clues to discover the actual nature of dark energy, whether it is
$\Lambda$CDM or another one. The analyses done in this work can be improved
with future more numerous and accurate data. Even though the $\Lambda$CDM
model was preferred, future higher redshift in GW data (due that GWTC-1 and
GWTC-2 do not exceed $z=1$) might change the preferences of GR over the siren
distance-luminosity distance coming from modified gravity models. This study
will be reported elsewhere.

Figure 7: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$ and $\Omega_{m}$ for the $\Lambda$CDM model. Left: Analyses
performed using Pantheon compilation (red color) and CC data (orange color).
Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 8: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in $\Lambda$CDM $h$, $\Omega_{M}$, $\Xi_{0}$ and $n$ using the
combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 9: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the CPL model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 10: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in CPL $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 11: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the BA model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 12: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in BA $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 13: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the JBP model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 14: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in JBP $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 15: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the LR model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 16: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in LR $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 17: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the LC model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 18: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in LC $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 19: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{m}$, $w_{0}$ and $w_{1}$ for the WP model. Left:
Analyses performed using Pantheon compilation (red color) and CC data (orange
color). Right: Analyses performed using Pantheon compilation, CC data and GW
catalogues (green color).

Figure 20: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters in WP $h$, $\Omega_{m}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ using
the combinations of samples: Pantheon+GW (black color), GW catalogues (green
color) and Pantheon+CC+GW (blue color). Left: case with $\delta_{1}$, middle:
case with $\delta_{2}$ and right: case with $\delta_{3}$.

Figure 21: 1-$\sigma$ and 2-$\sigma$ confidence contours for the cosmological
parameters $h$, $\Omega_{M}$, $w_{0}$, $w_{1}$, $\Xi_{0}$ and $n$ for the
seven models. The overlapping corresponds to the three delta functions plus
the GR case. From bottom left up to bottom: $\Lambda$CDM, CPL, BA, JBP, LR, LC
and WP, respectively.

Figure 22: Left: Siren distance-luminosity distance ratio for the $\Lambda$CDM model and the six dark energy parameterisations with the $\delta_{3}$ function with the $\Xi_{0}$ and $n$ best fits. The GR case is also plotted to compare with the other models. Right: $\delta_{3}$ function for the $\Lambda$CDM model and the six dark energy parameterisations with the $\delta_{3}$ function with the $\Xi_{0}$ and $n$ best fits. The GR case is also plotted to compare with the other models. Model | Database | $h$ | $\Omega_{M}$ | $w_{0}$ | $w_{1}$
---|---|---|---|---|---
$\Lambda CDM$ | CC | $0.7152\pm 0.0099$ | $0.248\pm 0.014$ | - | -
Pantheon | $0.7284\pm 0.0023$ | $0.285\pm 0.013$ | - | -
Pantheon+CC | $0.7330\pm 0.0020$ | $0.2447\pm 0.0073$ | - | -
CPL | CC | $0.722\pm 0.017$ | $0.186^{+0.071}_{-0.032}$ | $-1.05\pm 0.12$ | $0.77^{+0.54}_{-0.21}$
Pantheon | $0.7335\pm 0.0042$ | $0.306^{+0.13}_{-0.044}$ | $-1.18\pm 0.17$ | $-0.05^{+1.9}_{-0.61}$
Pantheon+CC | $0.7337\pm 0.0031$ | $0.168^{+0.078}_{-0.031}$ | $-0.950^{+0.053}_{-0.074}$ | $0.61^{+0.46}_{-0.17}$
BA | CC | $0.722\pm 0.016$ | $0.184^{+0.073}_{-0.031}$ | $-1.01^{+0.10}_{-0.12}$ | $0.380^{+0.25}_{-0.090}$
Pantheon | $0.7329\pm 0.0041$ | $0.336^{+0.11}_{-0.034}$ | $-1.20^{+0.15}_{-0.19}$ | $-0.40^{+1.3}_{-0.43}$
Pantheon+CC | $0.7330\pm 0.0029$ | $0.159^{+0.085}_{-0.038}$ | $-0.915^{+0.057}_{-0.090}$ | $0.326^{+0.20}_{-0.061}$
JBP | CC | $0.732^{+0.017}_{-0.015}$ | $0.156^{+0.081}_{-0.038}$ | $-1.37^{+0.14}_{-0.16}$ | $3.51^{+1.3}_{-0.42}$
Pantheon | $0.7344\pm 0.0045$ | $0.299^{+0.12}_{-0.051}$ | $-1.21^{+0.18}_{-0.16}$ | $0.49^{+2.5}_{-0.97}$
Pantheon+CC | $0.7347\pm 0.0039$ | $0.217^{+0.024}_{-0.018}$ | $-1.041\pm 0.079$ | $0.86\pm 0.73$
LR | CC | $0.718\pm 0.016$ | $0.226^{+0.046}_{-0.020}$ | $-0.992^{+0.094}_{-0.11}$ | $-0.09^{+0.12}_{-0.31}$
Pantheon | $0.7328\pm 0.0040$ | $0.324^{+0.12}_{-0.039}$ | $-1.18^{+0.16}_{-0.19}$ | $0.37^{+0.43}_{-1.5}$
Pantheon+CC | $0.7314\pm 0.0028$ | $0.211^{+0.044}_{-0.020}$ | $-0.936^{+0.043}_{-0.058}$ | $-0.096^{+0.095}_{-0.24}$
LC | CC | ${}^{\prime}0.723\pm 0.016$ | $0.188^{+0.070}_{-0.028}$ | $-1.05\pm 0.12$ | $-0.80^{+0.19}_{-0.16}$
Pantheon | $0.7336\pm 0.0041$ | $0.304^{+0.13}_{-0.047}$ | $-1.18\pm 0.17$ | $-1.18^{+0.75}_{-0.31}$
Pantheon+CC | $0.7337^{+0.0032}_{-0.0029}$ | $0.170^{+0.078}_{-0.031}$ | $-0.952^{+0.050}_{-0.074}$ | $-0.75^{+0.19}_{-0.15}$
WP | Chronometers | $0.729\pm 0.017$ | $0.105^{+0.044}_{-0.095}$ | $-1.13^{+0.14}_{-0.12}$ | $1.07^{+0.48}_{-0.36}$
Pantheon | $0.7352\pm 0.0043$ | $0.22^{+0.16}_{-0.14}$ | $-1.14^{+0.19}_{-0.15}$ | $0.66^{+0.77}_{-0.54}$
Pantheon+CC | $0.7354\pm 0.0035$ | $0.135^{+0.090}_{-0.050}$ | $-0.959\pm 0.054$ | $0.57^{+0.34}_{-0.31}$
Table 2: Best fit values for the four cosmological parameters with their
one-$\sigma$ uncertainties. The best fit values are reported for the standard
model $\Lambda$CDM and the six dynamical dark energy parameterisations. The
second column indicates the database used.

Model | Database (using $\delta_{1}$) | $h$ | $\Omega_{M}$ | $w_{0}$ | $w_{1}$ | $\Xi_{0}$ | $n$
---|---|---|---|---|---|---|---
$\Lambda CDM$ | GW | $0.690\pm 0.073$ | $0.26\pm 0.14$ | - | - | $1.00\pm 0.40$ | $1.93^{+0.77}_{-1.9}$
Pantheon+GW | $0.7283\pm 0.0023$ | $0.286\pm 0.013$ | - | - | $1.19^{+0.24}_{-0.30}$ | $1.8^{+1.2}_{-1.8}$
Pantheon+CC+GW | $0.7331\pm 0.0019$ | $0.2447\pm 0.0072$ | - | - | $1.17^{+0.25}_{-0.28}$ | $1.7^{+1.2}_{-1.7}$
CPL | GW | $0.698\pm 0.078$ | $0.286^{+0.21}_{-0.079}$ | $-1.4^{+2.2}_{-1.7}$ | $-0.1^{+2.9}_{-4.6}$ | $1.04^{+0.48}_{-0.57}$ | $2.2^{+1.5}_{-2.1}$
Pantheon+GW | $0.7335\pm 0.0041$ | $0.301^{+0.13}_{-0.046}$ | $-1.18\pm 0.17$ | $0.08^{+1.7}_{-0.62}$ | $1.18^{+0.24}_{-0.30}$ | $1.8^{+1.1}_{-1.8}$
Pantheon+CC+GW | $0.7337\pm 0.0031$ | $0.169^{+0.077}_{-0.030}$ | $-0.951^{+0.054}_{-0.074}$ | $0.61^{+0.46}_{-0.18}$ | $1.16\pm 0.34$ | $1.75^{+0.98}_{-1.8}$
BA | GW | $0.698\pm 0.078$ | $0.281^{+0.21}_{-0.093}$ | $-1.4^{+2.3}_{-1.6}$ | $-0.2^{+2.1}_{-4.4}$ | $1.06^{+0.51}_{-0.58}$ | $2.2^{+1.1}_{-2.2}$
Pantheon+GW | $0.7329\pm 0.0040$ | $0.332^{+0.11}_{-0.033}$ | $-1.20^{+0.15}_{-0.19}$ | $-0.34^{+1.3}_{-0.42}$ | $1.18^{+0.25}_{-0.29}$ | $1.8^{+1.6}_{-1.8}$
Pantheon+CC+GW | $0.7331\pm 0.0029$ | $0.157^{+0.088}_{-0.037}$ | $-0.913^{+0.056}_{-0.092}$ | $0.324^{+0.21}_{-0.056}$ | $1.17^{+0.24}_{-0.29}$ | $1.8^{+1.0}_{-1.8}$
JBP | GW | $0.698\pm 0.077$ | $0.28^{+0.21}_{-0.14}$ | $-1.4^{+2.2}_{-1.6}$ | $-0.1^{+3.8}_{-4.5}$ | $1.06^{+0.49}_{-0.56}$ | $2.17^{+0.87}_{-2.1}$
Pantheon+GW | $0.7344\pm 0.0044$ | $0.299^{+0.12}_{-0.044}$ | $-1.21^{+0.17}_{-0.15}$ | $0.6^{+2.3}_{-1.0}$ | $1.18^{+0.24}_{-0.28}$ | $1.8^{+1.1}_{-1.8}$
Pantheon+CC+GW | $0.7348\pm 0.0040$ | $0.216^{+0.024}_{-0.018}$ | $-1.043^{+0.074}_{-0.083}$ | $0.89\pm 0.74$ | $1.17^{+0.24}_{-0.29}$ | $1.77^{+0.69}_{-1.8}$
LR | GW | $0.698^{+0.085}_{-0.075}$ | $0.282^{+0.21}_{-0.092}$ | $-1.4^{+2.3}_{-1.6}$ | $0.2\pm 2.9$ | $1.07\pm 0.49$ | $2.2^{+1.2}_{-2.1}$
Pantheon+GW | $0.7328\pm 0.0040$ | $0.324^{+0.12}_{-0.034}$ | $-1.18^{+0.15}_{-0.20}$ | $0.33^{+0.48}_{-1.4}$ | $1.17^{+0.25}_{-0.30}$ | $1.75^{+0.69}_{-1.8}$
Pantheon+CC+GW | $0.7314\pm 0.0028$ | $0.211^{+0.043}_{-0.020}$ | $-0.937^{+0.044}_{-0.056}$ | $-0.098^{+0.095}_{-0.24}$ | $1.16^{+0.25}_{-0.28}$ | $1.77^{+0.94}_{-1.8}$
LC | GW | $0.691^{+0.098}_{-0.077}$ | $0.28^{+0.21}_{-0.13}$ | $-0.9^{+1.7}_{-3.3}$ | $-0.8^{+1.8}_{-3.4}$ | $1.12\pm 0.48$ | $2.18^{+0.88}_{-2.1}$
Pantheon+GW | $0.7333\pm 0.0042$ | $0.309^{+0.13}_{-0.044}$ | $-1.18\pm 0.17$ | $-1.21^{+0.77}_{-0.32}$ | $1.17^{+0.25}_{-0.29}$ | $1.8^{+1.4}_{-1.8}$
Pantheon+CC+GW | $0.7337\pm 0.0031$ | $0.172^{+0.075}_{-0.029}$ | $-0.953^{+0.051}_{-0.069}$ | $-0.75^{+0.18}_{-0.15}$ | $1.17^{+0.26}_{-0.29}$ | $1.74^{+0.61}_{-1.8}$
WP | GW | $0.695\pm 0.079$ | $0.27^{+0.22}_{-0.12}$ | $-1.1\pm 2.1$ | $2.3^{+2.6}_{-1.7}$ | $1.14^{+0.68}_{-0.44}$ | $2.2^{+1.6}_{-2.2}$
Pantheon+GW | $0.7351\pm 0.0044$ | $0.22\pm 0.12$ | $-1.14^{+0.19}_{-0.15}$ | $0.67^{+0.76}_{-0.58}$ | $1.17\pm 0.33$ | $1.8^{+1.0}_{-1.8}$
Pantheon+CC+GW | $0.7353\pm 0.0034$ | $0.134^{+0.086}_{-0.056}$ | $-0.958\pm 0.054$ | $0.58^{+0.36}_{-0.28}$ | $1.17^{+0.25}_{-0.29}$ | $1.8^{+1.1}_{-1.8}$

Table 3: Best fit values for the seven cosmological models discussed in Sec.IV
at 1-$\sigma$ uncertainties. In this case we consider the GW catalogues using
$\delta_{1}$ function (23).

Model | Database (using $\delta_{2}$) | $h$ | $\Omega_{M}$ | $w_{0}$ | $w_{1}$ | $\Xi_{0}$ | $n$
---|---|---|---|---|---|---|---
$\Lambda$CDM | GW | $0.686\pm 0.067$ | $0.26^{+0.22}_{-0.15}$ | - | - | $0.95^{+0.44}_{-0.56}$ | $2.72^{+2.2}_{-0.92}$
Pantheon+GW | $0.7283\pm 0.0023$ | $0.286\pm 0.013$ | - | - | $1.16^{+0.49}_{-0.35}$ | $2.92^{+2.0}_{-0.81}$
Pantheon+CC+GW | $0.7330\pm 0.0019$ | $0.2446\pm 0.0072$ | - | - | $1.14^{+0.55}_{-0.35}$ | $2.93^{+2.0}_{-0.81}$
CPL | GW | $0.696^{+0.095}_{-0.077}$ | $0.30^{+0.19}_{-0.11}$ | $-1.6^{+2.0}_{-1.4}$ | $-0.1\pm 2.9$ | $1.01\pm 0.52$ | $2.6\pm 1.4$
Pantheon+GW | $0.7334\pm 0.0042$ | $0.307^{+0.13}_{-0.043}$ | $-1.18\pm 0.17$ | $-0.01^{+1.8}_{-0.71}$ | $1.17^{+0.50}_{-0.36}$ | $2.91^{+2.0}_{-0.83}$
Pantheon+CC+GW | $0.7337\pm 0.0031$ | $0.167^{+0.079}_{-0.033}$ | $-0.949^{+0.053}_{-0.075}$ | $0.62^{+0.45}_{-0.16}$ | $1.13^{+0.51}_{-0.36}$ | $2.96^{+1.9}_{-0.88}$
BA | GW | $0.695^{+0.096}_{-0.073}$ | $0.297^{+0.20}_{-0.076}$ | $-1.6^{+2.0}_{-1.4}$ | $-0.1^{+2.6}_{-4.4}$ | $0.99^{+0.52}_{-0.63}$ | $2.6\pm 1.4$
Pantheon+GW | $0.7327\pm 0.0041$ | $0.331^{+0.12}_{-0.033}$ | $-1.19^{+0.15}_{-0.19}$ | $-0.36^{+1.3}_{-0.44}$ | $1.16^{+0.51}_{-0.35}$ | $2.89^{+2.0}_{-0.85}$
Pantheon+CC+GW | $0.7331\pm 0.0028$ | $0.157^{+0.086}_{-0.037}$ | $-0.914^{+0.052}_{-0.091}$ | $0.332^{+0.20}_{-0.062}$ | $1.14^{+0.48}_{-0.36}$ | $2.93^{+2.0}_{-0.79}$
JBP | GW | $0.696^{+0.093}_{-0.072}$ | $0.299^{+0.19}_{-0.078}$ | $-1.6^{+2.0}_{-1.3}$ | $0.0\pm 2.9$ | $1.01\pm 0.51$ | $2.6\pm 1.5$
Pantheon+GW | $0.7344\pm 0.0044$ | $0.296^{+0.12}_{-0.043}$ | $-1.21\pm 0.16$ | $0.6^{+2.3}_{-1.1}$ | $1.17^{+0.50}_{-0.36}$ | $2.91^{+2.0}_{-0.83}$
Pantheon+CC+GW | $0.7349\pm 0.0040$ | $0.216^{+0.024}_{-0.018}$ | $-1.045\pm 0.080$ | $0.89\pm 0.74$ | $1.14^{+0.49}_{-0.38}$ | $2.9^{+2.0}_{-1.1}$
LR | GW | $0.693^{+0.094}_{-0.077}$ | $0.298^{+0.19}_{-0.074}$ | $-1.6^{+2.0}_{-1.4}$ | $0.1\pm 2.9$ | $0.997\pm 0.53$ | $2.6\pm 1.5$
Pantheon+GW | $0.7327\pm 0.0040$ | $0.323^{+0.12}_{-0.037}$ | $-1.18^{+0.16}_{-0.20}$ | $0.33^{+0.47}_{-1.4}$ | $1.16^{+0.50}_{-0.35}$ | $2.91^{+2.0}_{-0.81}$
Pantheon+CC+GW | $0.7313\pm 0.0028$ | $0.210^{+0.045}_{-0.019}$ | $-0.935^{+0.044}_{-0.058}$ | $-0.101^{+0.094}_{-0.23}$ | $1.15^{+0.50}_{-0.38}$ | $2.9^{+2.0}_{-1.0}$
LC | GW | $0.687^{+0.10}_{-0.076}$ | $0.294^{+0.20}_{-0.092}$ | $-1.4^{+1.6}_{-2.9}$ | $-0.6^{+2.2}_{-3.8}$ | $1.06^{+0.64}_{-0.57}$ | $2.7^{+2.2}_{-1.4}$
Pantheon+GW | $0.7334\pm 0.0042$ | $0.310^{+0.13}_{-0.044}$ | $-1.18\pm 0.17$ | $-1.21^{+0.77}_{-0.30}$ | $1.15^{+0.51}_{-0.35}$ | $2.91^{+2.0}_{-0.79}$
Pantheon+CC+GW | $0.7337\pm 0.0031$ | $0.168^{+0.079}_{-0.028}$ | $-0.950^{+0.050}_{-0.075}$ | $-0.75\pm 0.15$ | $1.13^{+0.50}_{-0.39}$ | $2.93^{+2.0}_{-0.83}$
WP | GW | $0.687^{+0.093}_{-0.078}$ | $0.280^{+0.20}_{-0.088}$ | $-1.6\pm 1.8$ | $2.3^{+2.6}_{-1.1}$ | $1.06^{+0.68}_{-0.46}$ | $2.6\pm 1.5$
Pantheon+GW | $0.7350\pm 0.0043$ | $0.22^{+0.17}_{-0.14}$ | $-1.14^{+0.19}_{-0.15}$ | $0.66^{+0.77}_{-0.57}$ | $1.17^{+0.49}_{-0.37}$ | $2.95^{+1.9}_{-0.78}$
Pantheon+CC+GW | $0.7352\pm 0.0035$ | $0.135^{+0.087}_{-0.055}$ | $-0.957\pm 0.055$ | $0.57^{+0.36}_{-0.29}$ | $1.13^{+0.51}_{-0.39}$ | $2.96^{+2.0}_{-0.78}$

Table 4: Best fit values for the seven cosmological models discussed in Sec.IV
at 1-$\sigma$ uncertainties. In this case we consider the GW catalogues using
$\delta_{2}$ function (24).

Model | Database (using $\delta_{3}$) | $h$ | $\Omega_{M}$ | $w_{0}$ | $w_{1}$ | $\Xi_{0}$ | $n$
---|---|---|---|---|---|---|---
$\Lambda$CDM | GW | $0.694^{+0.054}_{-0.061}$ | $0.26^{+0.22}_{-0.12}$ | - | - | $1.04\pm 0.47$ | $1.81^{+0.59}_{-1.8}$
Pantheon+GW | $0.7283\pm 0.0023$ | $0.286\pm 0.013$ | - | - | $1.19^{+0.43}_{-0.37}$ | $1.85^{+0.61}_{-1.8}$
Pantheon+CC+GW | $0.7330\pm 0.0019$ | $0.2446\pm 0.0071$ | - | - | $1.17^{+0.47}_{-0.35}$ | $1.8^{+1.1}_{-1.8}$
CPL | GW | $0.695\pm 0.076$ | $0.295^{+0.20}_{-0.080}$ | $-1.5^{+2.0}_{-1.6}$ | $-0.1\pm 2.9$ | $1.07\pm 0.50$ | $2.2^{+2.5}_{-2.2}$
Pantheon+GW | $0.7335\pm 0.0042$ | $0.304^{+0.13}_{-0.045}$ | $-1.18\pm 0.17$ | $0.03^{+1.8}_{-0.68}$ | $1.19^{+0.49}_{-0.34}$ | $1.78^{+0.59}_{-1.7}$
Pantheon+CC+GW | $0.7336\pm 0.0031$ | $0.170^{+0.075}_{-0.029}$ | $-0.951^{+0.050}_{-0.072}$ | $0.60^{+0.46}_{-0.18}$ | $1.18^{+0.45}_{-0.38}$ | $1.79^{+0.77}_{-1.8}$
BA | GW | $0.696\pm 0.077$ | $0.29^{+0.20}_{-0.11}$ | $-1.5^{+2.1}_{-1.7}$ | $-0.1^{+2.7}_{-4.5}$ | $1.06^{+0.80}_{-0.58}$ | $2.1^{+1.3}_{-2.1}$
Pantheon+GW | $0.7328\pm 0.0039$ | $0.327^{+0.12}_{-0.036}$ | $-1.19^{+0.15}_{-0.20}$ | $-0.33^{+1.2}_{-0.43}$ | $1.19^{+0.47}_{-0.35}$ | $1.81^{+0.57}_{-1.8}$
Pantheon+CC+GW | $0.7330\pm 0.0029$ | $0.160^{+0.083}_{-0.036}$ | $-0.916^{+0.053}_{-0.087}$ | $0.323^{+0.20}_{-0.075}$ | $1.17^{+0.44}_{-0.39}$ | $1.78^{+0.57}_{-1.8}$
JBP | GW | $0.697\pm 0.076$ | $0.29^{+0.20}_{-0.12}$ | $-1.4^{+2.1}_{-1.6}$ | $-0.2\pm 2.9$ | $1.08\pm 0.51$ | $2.2\pm 1.4$
Pantheon+GW | $0.7344\pm 0.0044$ | $0.295^{+0.12}_{-0.047}$ | $-1.20\pm 0.16$ | $0.60^{+2.4}_{-0.98}$ | $1.19^{+0.45}_{-0.38}$ | $1.82^{+0.57}_{-1.8}$
Pantheon+CC+GW | $0.7347\pm 0.0040$ | $0.216^{+0.024}_{-0.017}$ | $-1.040^{+0.074}_{-0.085}$ | $0.87\pm 0.74$ | $1.17^{+0.48}_{-0.34}$ | $1.77^{+0.56}_{-1.8}$
LR | GW | $0.697\pm 0.078$ | $0.294^{+0.20}_{-0.074}$ | $-1.5^{+2.1}_{-1.6}$ | $0.1\pm 2.9$ | $1.07\pm 0.50$ | $2.16^{+0.95}_{-2.1}$
Pantheon+GW | $0.7328\pm 0.0040$ | $0.326^{+0.12}_{-0.037}$ | $-1.19^{+0.15}_{-0.19}$ | $0.36^{+0.49}_{-1.4}$ | $1.18^{+0.45}_{-0.37}$ | $1.82^{+0.57}_{-1.8}$
Pantheon+CC+GW | $0.7313\pm 0.0028$ | $0.210^{+0.044}_{-0.019}$ | $-0.935^{+0.040}_{-0.060}$ | $-0.102^{+0.092}_{-0.23}$ | $1.17^{+0.46}_{-0.37}$ | $1.74^{+0.52}_{-1.7}$
LC | GW | $0.686\pm 0.084$ | $0.294^{+0.20}_{-0.078}$ | $-1.2^{+1.6}_{-3.0}$ | $-0.6^{+1.9}_{-3.7}$ | $1.13\pm 0.50$ | $2.2^{+2.4}_{-2.1}$
Pantheon+GW | $0.7336\pm 0.0041$ | $0.305^{+0.13}_{-0.048}$ | $-1.18\pm 0.17$ | $-1.19^{+0.77}_{-0.29}$ | $1.18^{+0.47}_{-0.33}$ | $1.80^{+0.61}_{-1.8}$
Pantheon+CC+GW | $0.7336\pm 0.0030$ | $0.168^{+0.079}_{-0.032}$ | $-0.949^{+0.052}_{-0.075}$ | $-0.75^{+0.19}_{-0.15}$ | $1.18^{+0.44}_{-0.38}$ | $1.80^{+0.61}_{-1.8}$
WP | GW | $0.687\pm 0.079$ | $0.28^{+0.21}_{-0.11}$ | $-1.3^{+1.9}_{-2.3}$ | $2.3^{+2.5}_{-1.0}$ | $1.14^{+0.69}_{-0.43}$ | $2.2^{+1.8}_{-2.1}$
Pantheon+GW | $0.7352\pm 0.0044$ | $0.23^{+0.17}_{-0.10}$ | $-1.14^{+0.19}_{-0.15}$ | $0.66^{+0.75}_{-0.58}$ | $1.19^{+0.43}_{-0.38}$ | $1.85^{+0.71}_{-1.8}$
Pantheon+CC+GW | $0.7353\pm 0.0035$ | $0.135^{+0.090}_{-0.053}$ | $-0.958\pm 0.054$ | $0.57^{+0.37}_{-0.30}$ | $1.17^{+0.47}_{-0.37}$ | $1.80^{+0.59}_{-1.8}$

Table 5: Best fit values for the seven cosmological models discussed in Sec.IV at 1-$\sigma$ uncertainties. In this case we consider the GW catalogues using $\delta_{3}$ function (25). Model | $\chi_{\text{min}}^{2}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$BIC
---|---|---|---|---|---
$\Lambda$CDM | 1038.15 | 1046.15 | 1066.34 | - | -
CPL | 1035.24 | 1047.24 | 1077.52 | 1.08 | 11.17
BA | 1034.28 | 1046.28 | 1076.57 | 0.12 | 10.22
JBP | 1035.69 | 1047.69 | 1077.97 | 1.53 | 11.62
LR | 1036.42 | 1048.42 | 1078.70 | 2.26 | 12.35
LC | 1035.24 | 1047.24 | 1077.52 | 1.08 | 11.17
WP | 1035.27 | 1047.27 | 1077.55 | 1.11 | 11.20
Table 6: AIC and BIC criteria for the six models dark energy against $\Lambda$CDM for $\delta_{1}$ (23). Model | $\chi_{\text{min}}^{2}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$BIC
---|---|---|---|---|---
$\Lambda$CDM | 1038.66 | 1046.66 | 1066.85 | - | -
CPL | 1035.72 | 1047.72 | 1078.01 | 1.06 | 11.15
BA | 1034.77 | 1046.77 | 1077.05 | 0.10 | 10.19
JBP | 1036.17 | 1048.17 | 1078.46 | 1.51 | 11.60
LR | 1036.91 | 1048.91 | 1079.20 | 2.25 | 12.34
LC | 1035.72 | 1047.72 | 1078.01 | 1.06 | 11.15
WP | 1035.75 | 1047.75 | 1078.04 | 1.08 | 11.18
Table 7: AIC and BIC criteria for the six models dark energy against the $\Lambda$CDM for $\delta_{2}$ (24). Model | $\chi_{\text{min}}^{2}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$ BIC
---|---|---|---|---|---
$\Lambda$CDM | 1038.41 | 1046.41 | 1066.60 | - | -
CPL | 1035.36 | 1047.36 | 1077.65 | 0.95 | 11.05
BA | 1034.55 | 1046.55 | 1076.84 | 0.14 | 10.24
JBP | 1035.96 | 1047.96 | 1078.24 | 1.54 | 11.64
LR | 1036.69 | 1048.69 | 1078.97 | 2.28 | 12.37
LC | 1035.50 | 1047.50 | 1077.79 | 1.09 | 11.19
WP | 1035.52 | 1047.52 | 1077.81 | 1.11 | 11.20
Table 8: AIC and BIC criteria for the six models dark energy against the $\Lambda$CDM for $\delta_{3}$ (25). Model | $\delta$ | $\chi^{2}$ | AIC | BIC |  $\Delta$AIC |  $\Delta$BIC
---|---|---|---|---|---|---
$\Lambda$CDM | GR | 0.29 | 4.29 | 8.11 | - | -
$\delta_{1}$ | 0.29 | 8.29 | 15.94 | 3.99 | 7.82
$\delta_{2}$ | 0.29 | 8.29 | 15.94 | 3.99 | 7.82
$\delta_{3}$ | 0.28 | 8.28 | 15.93 | 3.98 | 7.81
CPL | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.29 | 12.29 | 23.76 | 3.99 | 7.82
$\delta_{2}$ | 0.28 | 12.28 | 23.75 | 3.99 | 7.81
$\delta_{3}$ | 0.28 | 12.28 | 23.76 | 3.99 | 7.81
BA | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.27 | 12.27 | 23.75 | 3.98 | 7.81
$\delta_{2}$ | 0.28 | 12.28 | 23.75 | 3.99 | 7.81
$\delta_{3}$ | 0.24 | 12.24 | 23.71 | 3.94 | 7.77
JBP | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.29 | 12.29 | 23.76 | 3.99 | 7.82
$\delta_{2}$ | 0.29 | 12.28 | 23.75 | 3.99 | 7.81
$\delta_{3}$ | 0.24 | 12.24 | 23.71 | 3.95 | 7.77
LR | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.28 | 12.28 | 23.75 | 3.99 | 7.81
$\delta_{2}$ | 0.28 | 12.28 | 23.75 | 3.99 | 7.81
$\delta_{3}$ | 0.28 | 12.28 | 23.76 | 3.99 | 7.81
LC | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.28 | 12.28 | 23.75 | 3.98 | 7.81
$\delta_{2}$ | 0.28 | 12.28 | 23.75 | 3.98 | 7.81
$\delta_{3}$ | 0.27 | 12.27 | 23.74 | 3.98 | 7.807
WP | GR | 0.29 | 8.29 | 15.94 | - | -
$\delta_{1}$ | 0.29 | 12.29 | 23.76 | 3.99 | 7.82
$\delta_{2}$ | 0.28 | 12.28 | 23.75 | 3.98 | 7.81
$\delta_{3}$ | 0.24 | 12.24 | 23.71 | 3.94 | 7.77
Table 9: AIC and BIC criteria for the three delta functions $\delta_{1}$ (23),
$\delta_{2}$ (24), $\delta_{3}$ (25) and the six dark energy models and the
$\Lambda$CDM against the GR model where $D_{L}(z)=D_{S}(z)$.

###### Acknowledgements.

CE-R acknowledges the Royal Astronomical Society as FRAS 10147. CE-R and AN
are supported by PAPIIT Project IA100220. This article is also based upon work
from COST action CA18108, supported by COST (European Cooperation in Science
and Technology). Numerical computations for this research were done in the
Tochtli ICN-UNAM cluster.

## References

  * [1] Celia Escamilla-Rivera.  Status on bidimensional dark energy parameterizations using SNe Ia JLA and BAO datasets.  Galaxies, 4(3):8, 2016.
  * [2] Celia Escamilla-Rivera, Luciano Casarini, Julio C. Fabris, and Jailson S. Alcaniz.  Linear and non-linear perturbations in dark energy models.  JCAP, 11:010, 2016.
  * [3] Luisa G. Jaime, Mariana Jaber, and Celia Escamilla-Rivera.  New parametrized equation of state for dark energy surveys.  Phys. Rev. D, 98(8):083530, 2018.
  * [4] Celia Escamilla-Rivera, A. Hernández-Almada, Miguel A. García-Aspeitia, and V. Motta.  Linear perturbations spectra for dynamical dark energy inspired by modified gravity.  arXiv preprint arXiv:2005.13957.
  * [5] N. M. Jiménez Cruz and Celia Escamilla-Rivera.  Late-time and Big Bang Nucleosynthesis constraints for generic modified gravity surveys.  Eur. Phys. J. Plus, 136(1):51, 2021.
  * [6] Sudipta Das, Abdulla Al Mamon, and Manisha Banerjee.  A new parametrization of dark energy equation of state leading to double exponential potential.  Res. Astron. Astrophys., 18(11):131, 2018.
  * [7] Emilio Elizalde, Martiros Khurshudyan, and Shin’ichi Nojiri.  Cosmological singularities in interacting dark energy models with an $\omega(q)$ parametrization.  Int. J. Mod. Phys. D, 28(01):1950019, 2018.
  * [8] Mehdi Rezaei, Mohammad Malekjani, Spyros Basilakos, Ahmad Mehrabi, and David F. Mota.  Constraints to Dark Energy Using PADE Parameterizations.  Astrophys. J., 843(1):65, 2017.
  * [9] Celia Escamilla-Rivera and Salvatore Capozziello.  Unveiling cosmography from the dark energy equation of state.  Int. J. Mod. Phys. D, 28(12):1950154, 2019.
  * [10] Cristian Zamora Munõz and Celia Escamilla-Rivera.  Inverse Cosmography: testing the effectiveness of cosmographic polynomials using machine learning.  JCAP, 12:007, 2020.
  * [11] Eleonora Di Valentino, Luis A Anchordoqui, Yacine Ali-Haimoud, Luca Amendola, Nikki Arendse, Marika Asgari, Mario Ballardini, Elia Battistelli, Micol Benetti, Simon Birrer, et al.  Cosmology intertwined i: Perspectives for the next decade.  arXiv preprint arXiv:2008.11283, 2020.
  * [12] B. P. Abbott et al.  Multi-messenger Observations of a Binary Neutron Star Merger.  Astrophys. J. Lett., 848(2):L12, 2017.
  * [13] Jose María Ezquiaga and Miguel Zumalacárregui.  Dark Energy in light of Multi-Messenger Gravitational-Wave astronomy.  Front. Astron. Space Sci., 5:44, 2018.
  * [14] Lorenzo Bordin, Edmund J. Copeland, and Antonio Padilla.  Dark energy loopholes some time after GW170817.  JCAP, 11:063, 2020.
  * [15] B. P. Abbott et al.  GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs.  Phys. Rev. X, 9(3):031040, 2019.
  * [16] Michel Chevallier and David Polarski.  Accelerating universes with scaling dark matter.  Int. J. Mod. Phys. D, 10:213–224, 2001.
  * [17] Yun Wang.  Figure of merit for dark energy constraints from current observational data.  Physical Review D, 77(12):123525, 2008.
  * [18] David W Hogg.  Distance measures in cosmology.  arXiv preprint astro-ph/9905116, 1999.
  * [19] S. Weinberg and W. Steven.  Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity.  Wiley, 1972.
  * [20] Scott Dodelson and Fabian Schmidt.  Modern Cosmology.  Academic Press, 2020.
  * [21] Bernard Schutz.  A first course in general relativity.  Cambridge university press, 2009.
  * [22] Sanjeev Dhurandhar and Bangalore Sathyaprakash.  Cosmic sirens: discovery of gravitational waves and their impact on astrophysics and fundamental physics.  Current Science, 113(4):663–671, 2017.
  * [23] Thibault Damour.  Gravitational waves and binary black holes.  Ondes Gravitationnelles, Séminaire Poincaré, 22:1–51, 2016\.
  * [24] Benjamin P Abbott, Richard Abbott, TD Abbott, MR Abernathy, Fausto Acernese, Kendall Ackley, Carl Adams, Thomas Adams, Paolo Addesso, RX Adhikari, et al.  Observation of gravitational waves from a binary black hole merger.  Physical review letters, 116(6):061102, 2016.
  * [25] Ayan Mitra, Jurgen Mifsud, David F. Mota, and David Parkinson.  Cosmology with the Einstein Telescope: No Slip Gravity Model and Redshift Specifications.  arXiv preprint arXiv:2010.00189.
  * [26] Enis Belgacem, Gianluca Calcagni, Marco Crisostomi, Charles Dalang, Yves Dirian, Jose María Ezquiaga, Matteo Fasiello, Stefano Foffa, Alexander Ganz, Juan García-Bellido, et al.  Testing modified gravity at cosmological distances with lisa standard sirens.  Journal of Cosmology and Astroparticle Physics, 2019(07):024, 2019\.
  * [27] Steven Weinberg.  The Cosmological constant problems.  In 4th ISSDDM (DM 2000), pages 18–26, 2 2000.
  * [28] Timothy Clifton, Pedro G. Ferreira, Antonio Padilla, and Constantinos Skordis.  Modified Gravity and Cosmology.  Phys. Rept., 513:1–189, 2012.
  * [29] Eric V. Linder.  The Dynamics of Quintessence, The Quintessence of Dynamics.  Gen. Rel. Grav., 40:329–356, 2008.
  * [30] EM Barboza Jr and JS Alcaniz.  A parametric model for dark energy.  Physics Letters B, 666(5):415–419, 2008.
  * [31] HK Jassal, JS Bagla, and T Padmanabhan.  Wmap constraints on low redshift evolution of dark energy.  Monthly Notices of the Royal Astronomical Society: Letters, 356(1):L11–L16, 2005.
  * [32] Dragan Huterer and Michael S. Turner.  Probing dark energy: Methods and strategies.  Phys. Rev. D, 64:123527, Nov 2001.
  * [33] Christof Wetterich.  Phenomenological parameterization of quintessence.  Physics Letters B, 594(1-2):17–22, 2004.
  * [34] Daniel Moshe Scolnic, DO Jones, A Rest, YC Pan, R Chornock, RJ Foley, ME Huber, R Kessler, Gautham Narayan, AG Riess, et al.  The complete light-curve sample of spectroscopically confirmed sne ia from pan-starrs1 and cosmological constraints from the combined pantheon sample.  The Astrophysical Journal, 859(2):101, 2018.
  * [35] David W Hogg, Jo Bovy, and Dustin Lang.  Data analysis recipes: Fitting a model to data.  arXiv preprint arXiv:1008.4686, 2010.
  * [36] Juan Magaña, Mario H Amante, Miguel A Garcia-Aspeitia, and V Motta.  The Cardassian expansion revisited: constraints from updated Hubble parameter measurements and type Ia supernova data.  Monthly Notices of the Royal Astronomical Society, 476(1):1036–1049, 02 2018.
  * [37] BP Abbott, Richard Abbott, TD Abbott, S Abraham, F Acernese, K Ackley, C Adams, RX Adhikari, VB Adya, C Affeldt, et al.  Gwtc-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs.  Physical Review X, 9(3):031040, 2019.
  * [38] R Abbott, TD Abbott, S Abraham, F Acernese, K Ackley, A Adams, C Adams, RX Adhikari, VB Adya, C Affeldt, et al.  Gwtc-2: Compact binary coalescences observed by ligo and virgo during the first half of the third observing run.  arXiv preprint arXiv:2010.14527, 2020.
  * [39] Daniel Foreman-Mackey, David W Hogg, Dustin Lang, and Jonathan Goodman.  emcee: the mcmc hammer.  Publications of the Astronomical Society of the Pacific, 125(925):306, 2013.
  * [40] Antony Lewis.  Getdist: a python package for analysing monte carlo samples.  arXiv preprint arXiv:1910.13970, 2019.
  * [41] Benjamin Audren, Julien Lesgourgues, Karim Benabed, and Simon Prunet.  Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code.  JCAP, 02:001, 2013.
  * [42] Thejs Brinckmann and Julien Lesgourgues.  MontePython 3: boosted MCMC sampler and other features.  Phys. Dark Univ., 24:100260, 2019.
  * [43] Andrew Gelman and Donald B. Rubin.  Inference from Iterative Simulation Using Multiple Sequences.  Statist. Sci., 7:457–472, 1992.
  * [44] Harold Jeffreys.  The Theory of Probability.  Oxford Classic Texts in the Physical Sciences. 1939.
  * [45] Alan Heavens, Yabebal Fantaye, Arrykrishna Mootoovaloo, Hans Eggers, Zafiirah Hosenie, Steve Kroon, and Elena Sellentin.  Marginal Likelihoods from Monte Carlo Markov Chains.  arXiv preprint 1704.03472.
  * [46] Andrew R. Liddle.  How many cosmological parameters?  Mon. Not. Roy. Astron. Soc., 351:L49–L53, 2004.
  * [47] Robert E. Kass and Adrian E. Raftery.  Bayes factors.  Journal of the American Statistical Association, 90(430):773–795, 1995.

##  Appendix A Derivation of the tensor for Gravitational Waves

We can compute the $h_{ij}$ (the spacial part of $h_{\mu\nu}$) tensor using
the spatial part of the wave equation for GWs [21]

$\left(\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}-\nabla^{2}\right)h_{\mu\nu}=\dfrac{16\pi G}{c^{4}}T_{\mu\nu}.$ (59)

To achieve this, it is easier to work with the Green’s function. The Green’s
function equation for this can be expressed as

$\left(\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}-\nabla^{2}\right)G(\mathbf{x},\mathbf{x}_{0},t,t_{0})=\delta(\mathbf{x}-\mathbf{x}_{0})\delta(t-t_{0}).$
(60)

By writing the Fourier transform of the Green’s function and the Dirac’s
deltas, and solving the equation the Green’s function is

$G(\mathbf{x},\mathbf{x}_{0},t,t_{0})=\dfrac{\delta\left(t-t_{0}-\dfrac{|\mathbf{x}-\mathbf{x}_{0}|}{c}\right)}{4\pi|\mathbf{x}-\mathbf{x}_{0}|},$
(61)

then, we can obtain the following solution (here we only focus on the spacial
part of the solution)

$h_{ij}^{TT}=\frac{16\pi
G}{c^{4}}\bigints_{V}d^{3}x_{0}\,dt_{0}\,T_{ij}(\mathbf{x}_{0},t_{0})\dfrac{\delta\left(t-t_{0}-\dfrac{|\mathbf{x}-\mathbf{x}_{0}|}{c}\right)}{4\pi|\mathbf{x}-\mathbf{x}_{0}|},$
(62)

by integrating over $t_{0}$

$h_{ij}^{TT}=\frac{4G}{c^{4}}\bigints_{V}d^{3}x_{0}\frac{T_{ij}\left(t-\dfrac{|\mathbf{x}-\mathbf{x}_{0}|}{c},\mathbf{x}_{0}\right)}{|\mathbf{x}-\mathbf{x}_{0}|}.$
(63)

If the system is far away then $|\mathbf{x}-\mathbf{x}_{0}|\approx D_{S}$,
where $D_{S}$ is the siren distance to the source (the distance associated
with the GW is the siren distance), then

$h_{ij}^{TT}=\frac{4G}{c^{4}D_{S}}\int
d^{3}x_{0}T_{ij}\left(t-\frac{D_{S}}{c},\mathbf{x}_{0}\right).$ (64)

Due the gauge condition $\partial_{\nu}T^{\mu\nu}=0$ [21], the expressions
$\partial T^{\mu 0}/\partial t=-\partial T^{\mu i}/\partial x^{i}$ and
$\partial T^{00}/\partial t=-\partial T^{0i}/\partial x^{i}$ hold. By using
these expressions and the conditions under the source is far away (the surface
integrals that include the stress-energy tensor are equal to zero) and slow-
motion $\left(\partial x^{i}/\partial t\approx 0\right)$ holds then the
following expressions can be derived

$\displaystyle\int d^{3}xT_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\dfrac{d}{dt}\int
d^{3}x(T_{i0}x_{j}+T_{j0}x_{i}),\quad\quad\quad$ (65) $\displaystyle\int
d^{3}x(T_{i0}x_{j}+T_{j0}x_{i})$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\int d^{3}xT_{00}x_{i}x_{j}.$ (66)

These are obtained by applying integration by parts in the right side above
expressions and consider zero over some terms with the conditions explained in
the previous paragraph. Then by defining the quadrupole moment tensor by
$I_{ij}=\displaystyle\int d^{3}x\rho(\mathbf{x})x_{i}x_{j}$ the final
expression is

$h_{ij}=\frac{2G}{c^{4}D_{S}}\left[\frac{d^{2}I_{ij}\left(t-\dfrac{R}{c}\right)}{dt^{2}}\right]^{TT},$
(67)

where $TT$ denotes that is the transverse and traceless part of $I_{ij}$.