# Testing $\Lambda$-Free $f(Q)$ Cosmology

José Ferreira<EMAIL_ADDRESS>Instituto de Astrofísica e
Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa, Campo
Grande, Edificio C8, P-1749-016, Lisboa, Portugal  Departamento de Física,
Faculdade de Ciências da Universidade de Lisboa, Campo Grande, Edificio C8,
P-1749-016, Lisboa, Portugal  Tiago Barreiro<EMAIL_ADDRESS>Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da
Universidade de Lisboa, Campo Grande, Edificio C8, P-1749-016, Lisboa,
Portugal  ECEO, Universidade Lusófona de Humanidades e Tecnologias, Campo
Grande, 376, 1749-024 Lisboa, Portugal  José P. Mimoso
<EMAIL_ADDRESS>Instituto de Astrofísica e Ciências do Espaço,
Faculdade de Ciências da Universidade de Lisboa, Campo Grande, Edificio C8,
P-1749-016, Lisboa, Portugal  Departamento de Física, Faculdade de Ciências da
Universidade de Lisboa, Campo Grande, Edificio C8, P-1749-016, Lisboa,
Portugal  Nelson J. Nunes<EMAIL_ADDRESS>Instituto de
Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de
Lisboa, Campo Grande, Edificio C8, P-1749-016, Lisboa, Portugal  Departamento
de Física, Faculdade de Ciências da Universidade de Lisboa, Campo Grande,
Edificio C8, P-1749-016, Lisboa, Portugal

###### Abstract

We study a model of Symmetric Teleparallel gravity that is able to account for
the current accelerated expansion of the universe without the need for dark
energy component. We investigate this model by making use of dynamical system
analysis techniques to identify the regions of the parameter space with viable
cosmologies and constrain it using type Ia supernova (SnIa), cosmic microwave
background (CMB) data and make forecasts using standard siren (SS) events. We
conclude that this model is disfavored with respect to $\Lambda$CDM and
forthcoming standard siren events can be decisive in testing the viability of
the model.

##  1 Introduction

The standard model of cosmology, referred to as $\Lambda$CDM, provides us with
the most accurate description of the Universe at large scales. However,
despite its successes, it is not devoid of difficulties [1]. One of the most
debated of its shortcomings is the origin of the cosmological constant,
introduced to explain the current accelerated expansion of the Universe. Also,
tensions in the value of the expansion rate obtained from low redshift against
high redshift data, motivate the search for alternative descriptions.

Recently, it has been shown that two gravitational theories equivalent to
General Relativity (GR) exist, one based solely on non-metricity and the other
based solely on torsion. This equivalence is referred as the geometrical
trinity of gravity [2]. Generalization of these formulations have been
considered as $f(T)$, the teleparallel gravity, or $f(Q)$, the symmetric
teleparallel gravity [3, 4, 5]. These theories are particularly interesting as
even if they only produce slight deviations from $\Lambda$CDM at the
background level [6], they may however, give rise to signatures in the
evolution of the density perturbations [7] and in the propagation of
gravitational waves [8], potentially being able to solve the challenges that
$\Lambda$CDM faces [9, 10, 11, 12, 13, 14, 15, 16].

With the development of new techniques in gravitational wave (GW) astronomy,
cosmologists will soon have access to a new and exciting source of data:
standard siren (SS) events, characterized by the direct measurement of both
the luminosity distance and the redshift of merging compact objects.
Unfortunately, up to the current date, only GW170817 [17] has been confirmed
as a SS event and a suggested electromagnetic (EM) counterpart to GW190521
[18] was proposed in [19].

In this work, we focus on a model of $f(Q)$ gravity able to account for the
late time accelerated expansion of the Universe without adding a cosmological
constant. We develop a dynamical system analysis in order to look for regions
in parameter space that yield viable cosmologies, and constrain this model by
making use of both high and low redshift observables, namely, type Ia
supernova (SnIa), the cosmic microwave background (CMB) and in addition, make
forecasts using mock catalogs of SS events made for Laser Interferometer Space
Antenna (LISA) [20, 21] and the Einstein Telescope (ET). An online repository
complementary to this analysis is publicly available at [22].

This work is organized as follows: In section 2 we provide a brief
introduction to the formalism behind $f(Q)$ cosmology and introduce the
specific model we will be working with. In section 3, we perform a dynamical
system analysis applied to the model to identify the regions of the phase
space leading to late time acceleration. In section 4 we introduce the
datasets we use. In section 5 we show the constraints given by the datasets
introduced before. Finally, in section 6, we make an overview of this work.

##  2 Cosmology in $f(Q)$

In this work, we consider a theory of modified gravity based on an arbitrary
function of the non-metricity scalar. The action for such a theory of gravity
reads

$S=\int\sqrt{-g}\left[-\frac{c^{4}}{16\pi
G}f(Q)+\mathcal{L}_{m}\right]d^{4}x\,,$ (1)

where $Q$ is the non-metricity scalar that is defined as [23]

$Q=-\frac{1}{4}Q_{\alpha\beta\gamma}Q^{\alpha\beta\gamma}+\frac{1}{2}Q_{\alpha\beta\gamma}Q^{\gamma\beta\alpha}+\frac{1}{4}Q_{\alpha}Q^{\alpha}-\frac{1}{2}Q_{\alpha}\tilde{Q}^{\alpha}\,,$
(2)

with $Q_{\alpha\mu\nu}=\nabla_{\alpha}\,g_{\mu\nu}$ and the two independent
contractions of the non-metricity tensor are

$Q_{\mu}=Q_{\mu}{}^{\alpha}{}_{\alpha}\,,\hskip
28.45274pt\tilde{Q}^{\mu}=Q_{\alpha}{}^{\alpha\mu}\,.$ (3)

When varying the action in eq. 1 with respect to the metric, one obtains the
field equations [24]

$\displaystyle\frac{2}{\sqrt{-g}}\nabla_{\alpha}\left(\sqrt{-g}f_{Q}P^{\alpha\mu}{}_{\nu}\right)+\frac{1}{2}\delta^{\mu}_{\nu}f$
$\displaystyle+f_{Q}P^{\mu\alpha\beta}Q_{\nu\alpha\beta}=\frac{8\pi
G}{c^{4}}T^{\mu}{}_{\nu}\,,$ (4)

where $f_{Q}$ represents the partial derivative of $f(Q)$ with respect to $Q$.
$P^{\alpha\mu\nu}$, referred to as the non-metricity conjugate, is defined as

$P^{\alpha}{}_{\mu\nu}=-\frac{1}{2}L^{\alpha}{}_{\mu\nu}+\frac{1}{4}(Q^{\alpha}-\tilde{Q}^{\alpha})-\frac{1}{4}\delta^{\alpha}_{(\mu}Q_{\nu)}\,.$
(5)

where $L^{\alpha}{}_{\mu\nu}$ is the disformation tensor

$L^{\lambda}{}_{\mu\nu}=\frac{1}{2}g^{\lambda\beta}\left(-Q_{\mu\beta\nu}-Q_{\nu\beta\mu}+Q_{\beta\mu\nu}\right)\,.$
(6)

We use the coincident gauge to specify the affine connection [3, 25, 23] and
restrict our analysis to a flat, homogeneous and isotropic universe, described
by the FLRW metric

$ds^{2}=-c^{2}dt^{2}+a^{2}(t)(dx^{2}+dy^{2}+dz^{2})\,,$ (7)

where $a(t)$ is the scale factor of the universe, and $t$ is cosmic time $t$.
Under these assumptions, the non-metricity scalar defined in eq. 2 simplifies
to [24]

$Q=6H^{2}\,,$ (8)

where $H=\dot{a}/a$ is the Hubble function. Considering a universe permeated
by a perfect fluid composed of both matter and radiation, the modified first
Friedmann equation reads [24]

$6f_{Q}H^{2}-\frac{1}{2}f=8\pi G(\rho_{m}+\rho_{r})\,.$ (9)

As stated in [26], the covariant derivative of the stress-energy momentum
tensor is still zero, meaning that each component of the fluid obeys the
standard continuity equation

$\dot{\rho}+3\frac{\dot{a}}{a}\left(\rho+\frac{P}{c^{2}}\right)=0\,.$ (10)

where $P$ stands for pressure. This means that the scaling for each component
is the same as in $\Lambda$CDM and we can rewrite eq. 9 as

$\frac{2Qf_{Q}-f}{Q_{0}}=\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}}\,,$
(11)

where $\Omega_{i}=8\pi G\rho_{i,0}/3H_{0}^{2}$ and $E\equiv H/H_{0}$, with the
index $0$ denoting the value of that quantity today.

This modification of gravity also changes the equation of motion for the
propagation of gravitational waves [24]

$\bar{h}_{A}^{\prime\prime}+2\mathcal{H}(1+\delta(z))\bar{h}_{A}^{\prime}+k^{2}\bar{h}_{A}=0\,,$
(12)

where $A=\times,+$ represent the polarizations of the GW. The previous
equation differs from the case of $\Lambda$CDM by an additional friction
factor, $\delta(z)$, that for an $f(Q)$ model takes the form

$\delta(z)=\frac{d\ln{f_{Q}}}{2\mathcal{H}d\eta}\,.$ (13)

Following [8], this leads to a modification to the luminosity distance for a
GW event, such that

$d_{\text{GW}}(z)=\sqrt{\frac{f_{Q}^{(0)}}{f_{Q}}}\,\,d_{L}(z)\,,$ (14)

where $f_{Q}^{(0)}$ is the function $f_{Q}$ computed at the present day and
$d_{L}(z)$ is the standard luminosity distance function.

###  2.1 $f(Q)$ Cosmological Model

In this work we consider a model that differs from the standard interpretation
of gravity, i.e., GR, for which $f(Q)=Q$, by taking a multiplicative
exponential term of the form [27],

$f(Q)=Qe^{\lambda Q_{0}/Q}\,,$ (15)

where $\lambda$ is a constant, however not a free parameter. As we will see in
section 3, $\lambda$ is related to the abundance of matter today and has to be
positive for a viable universe.

By making use of eq. 8, we can see that the extra multiplicative term in this
model is an exponential that is inversely proportional to the Hubble function.
Given that the Hubble function increase with redshift, this model
exponentially approaches the case where $f(Q)=Q$, meaning that it has GR as a
limit case at early times.

Inserting the specific form of $f(Q)$ into eq. 11 we obtain the modified first
Friedmann equation for this model, which reads

$(E^{2}-2\lambda)e^{\lambda/E^{2}}=\Omega_{m}(1+z)^{3}+\Omega_{r}(1+z)^{4}\,,$
(16)

where $E(z)\equiv H(z)/H_{0}$.

As hinted before, this model falls back to GR in the limit of high redshifts.
For $z\gg 1$, eq. 16 is approximately equal to

$E^{2}\approx\Omega_{m}(1+z)^{3}+\Omega_{r}(1+z)^{4}+\lambda\,,$ (17)

which corresponds to a $\Lambda$CDM universe with $\Omega_{\Lambda}=\lambda$,
where its effect is negligible at these early times. When evaluated today
$(E^{2}\approx 1)$, eq. 16 reduces to

$E^{2}\approx
e^{-\lambda}\Omega_{m}(1+z)^{3}+e^{-\lambda}\Omega_{r}(1+z)^{4}+2\lambda\,,$
(18)

which corresponds once again to a $\Lambda$CDM universe, but now with
$\Omega_{\Lambda}=2\lambda$ and a correction in the relative abundance of
matter and radiation. We note that, in the absence of matter fields,
$\rho_{m}=\rho_{r}=0$, we obtain a de Sitter universe with $\lambda=1/2$.

If we insert the form of the function $f(Q)$ in eq. 14, then we obtain the
luminosity distance function for this specific model which now reads

$d_{\text{GW}}(z)=\sqrt{\frac{1-\lambda}{1-\lambda/E^{2}}}e^{\frac{\lambda}{2}(1-1/E^{2})}d_{L}(z)\,,$
(19)

This model was originally proposed in [27], where it was shown to be
statistically equivalent to $\Lambda$CDM for tests using low redshift
measurements. Additionally, it has been shown to pass Big Bang nucleosynthesis
constraints [28].

##  3 Dynamical System Analysis

In order to perform our dynamical systems analysis, we can start by re-writing
eq. 16 as

$1=\frac{e^{-\lambda/E^{2}}}{E^{2}}\left(\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}}\right)+\frac{2\lambda}{E^{2}}\,.$
(20)

Defining the quantities

$\displaystyle
x_{1}\equiv\frac{\Omega_{m}}{E^{2}e^{\lambda/E^{2}}a^{3}}\,,\hskip
7.11317ptx_{2}\equiv\frac{\Omega_{r}}{E^{2}e^{\lambda/E^{2}}a^{4}}\,,\hskip
7.11317ptx_{3}\equiv\frac{2\lambda}{E^{2}}\,,$ (21)

where $x_{1}$ is related to the evolution of the matter density in the
universe, $x_{2}$ with the radiation density and $x_{3}$ with the parameter
$\lambda$, they satisfy

$x_{1}+x_{2}+x_{3}=1\,.$ (22)

This equation shows that one of our coordinates in phase space is fully
determined by the value of the other two. As such, without loss of generality,
we choose to work with $x_{1}$ and $x_{2}$ alone.

Differentiating both $x_{1}$ and $x_{2}$ with respect to the number of
$e$-folds, $N=\ln{a}$, we obtain

$\displaystyle
x_{1}^{\prime}=-x_{1}\left((1+x_{1}+x_{2})\frac{E^{\prime}}{E}+3\right)\,,$
(23) $\displaystyle
x_{2}^{\prime}=-x_{2}\left((1+x_{1}+x_{2})\frac{E^{\prime}}{E}+4\right)\,,$
(24)

where the prime denotes the derivative with respect to $N$ and the ratio
$E^{\prime}/E$ can be expressed as a function of both $x_{1}$ and $x_{2}$ by

$\frac{E^{\prime}}{E}=\frac{-(3x_{1}+4x_{2})}{2-(x_{1}+x_{2})+(x_{1}+x_{2})^{2}}\,.$
(25)

###  3.1 Fixed Points

For this dynamical system, we find three distinct fixed points. Their
locations and corresponding stability is presented in table 1.

Point | Type | $x_{1}$ | $x_{2}$ | Stability | Eigenvalues
---|---|---|---|---|---
I |  $\lambda$ dominated | 0 | 0 | Stable | $(-4,-3)$
II | Matter dominated | 1 | 0 | Saddle | $(3,-1)$
III | Radiation dominated | 0 | 1 | Unstable | $(4,1)$
Table 1: Fixed points and corresponding stability for the dynamical system.

Fixed point I corresponds to a $\lambda$ dominated regime and is the only
stable fixed point in our system. Located at $x_{1},x_{2}=0$, corresponding to

$E^{2}=2\lambda\,,$ (26)

and therefore it only exists if $\lambda$ is positive.

By integrating the previous equation we obtain the evolution of the scale
factor as a function of time

$a=e^{\sqrt{2\lambda}H_{0}t}\,,$ (27)

meaning that we have an exponential expansion, corresponding to a point of
eternal inflation.

The fixed point II corresponds to a matter dominated regime located at
$x_{1}=1$ and $x_{2}=0$, implying that $\lambda=0$. It is a saddle node with
attracting trajectories along the direction of fixed point III. Setting
$\lambda=0$ in eq. 15, we have $f(Q)=Q$, corresponding to the GR limit, more
specifically a cold dark matter (CDM) universe, since this model does not
include a cosmological constant.

As for fixed point III, it corresponds to a radiation dominated regime and is
located at $x_{1}=0$ and $x_{2}=1$. It is an unstable node and like fixed
point II it is located in a region where $\lambda=0$, meaning that it is in
the GR limit with $f(Q)=Q$.

###  3.2 Trajectories

Solving the equations of motion eqs. 23 and 24 numerically, we obtain the
trajectories in phase space presented in fig. 1. From this figure, we can see
the existence of three disjoint regions related with the sign of $\lambda$:

  * •

$\lambda>0$: region colored in light blue, corresponding to the triangle
$x_{1}+x_{2}<1$. This is where the stable fixed point I is located, therefore
leading all trajectories inside this region towards a universe in eternal
inflation.

  * •

$\lambda=0$: region in light brown, corresponding to the line $x_{1}+x_{2}=1$
and a CDM universe. In this region the trajectories flow from the unstable
fixed point III, corresponding to a radiation dominated epoch, towards the
saddle fixed point II, corresponding to a matter dominated epoch.

  * •

$\lambda<0$: region in dark blue, corresponding to the plane $x_{1}+x_{2}>1$.
It has no fixed points and all trajectories diverge towards infinity, with $E$
asymptotically approaching zero.

Figure 1: Stream plot of the phase space for the dynamical system where fixed
point I (stable) is represented in green, fixed point II (saddle) in orange
and fixed point III (unstable) in red. The light blue marks the region where
$\lambda>0$, the light brown region corresponds to $\lambda=0$ and the dark
blue region to $\lambda<0$.

In order to agree with the observations of an expanding universe the value of
$\lambda$ must be positive. This is what was obtained in [27] where this model
was constrained using SnIa, baryonic acoustic oscillations, cosmic
chronometers, and redshift space distortions.

Qualitatively, the previous dynamical system is similar to a $\Lambda$CDM
universe where the value of the cosmological constant can be positive,
negative or zero [29]. However, in $\Lambda$CDM with a negative cosmological
constant, the universe eventually collapses into a Big Crunch.

###  3.3 Computing the value of $\lambda$

In order to compute the value of $\lambda$, we evaluate eq. 16 at the present
day, obtaining

$\left(\lambda-\frac{1}{2}\right)e^{\lambda-1/2}=-\frac{\Omega_{m}+\Omega_{r}}{2\sqrt{e}}\,.$
(28)

This equation only possesses a solution provided that
$\Omega_{m}+\Omega_{r}\leq 2e^{-1/2}$. This is not immediately ensured because
the sum of the densities is not necessarily equal to 1. The two possible
solutions for $\lambda$ are

$\displaystyle\lambda_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2}+W_{0}\left(-\frac{\Omega_{m}+\Omega_{r}}{2e^{1/2}}\right)\,,$
(29) $\displaystyle\lambda_{-1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}+W_{-1}\left(-\frac{\Omega_{m}+\Omega_{r}}{2e^{1/2}}\right)\,,$
(30)

where $W_{0}$ and $W_{-1}$ are the main and the $-1$ branch of the Lambert
function, respectively.

We discard the solution $\lambda_{-1}$, since it only provides negative values
for $\lambda$. Additionally, we see that $\lambda_{0}$ takes negative values
when $\Omega_{m}+\Omega_{r}>1$, hence, for $\lambda$ to be positive, the sum
of these densities must be smaller or equal to 1 as in $\Lambda$CDM.

##  4 Datasets

Throughout this work we make use of two different sources of observational
data, SnIa and CMB, and of mock catalogs for SS events generated for the LISA
space mission and ET.

###  4.1 Type Ia Supernovae

We use SnIa to provide low redshift constraints for this model. The dataset of
SnIa used is the Pantheon sample [30], with the corresponding data publicly
available in [31]. For performance reasons, the binned sample of this dataset
was used.

The relationship between the apparent magnitude and the luminosity distance is
given by

$m=M+5\log{\left(\frac{d_{L}(z)}{\rm Mpc}\right)}+25\,,$ (31)

where $M$ is the bolometric magnitude. In order to avoid degeneracies and
unwanted parameters, we followed the procedure presented in [32], where a
marginalization was performed on both $H_{0}$ and the bolometric magnitude.
The likelihood for the SnIa reads

$L=\exp\left[-\frac{1}{2}\left(A-\frac{B^{2}}{C}\right)\right]\,,$ (32)

where $A$, $B$ and $C$ are defined as

$\displaystyle
A=\sum_{i=1}^{n}\frac{\Delta^{2}(z_{i})}{\sigma^{2}(z_{i})}\,,B=\sum_{i=1}^{n}\frac{\Delta(z_{i})}{\sigma^{2}(z_{i})}\,,C=\sum_{i=1}^{n}\frac{1}{\sigma^{2}(z_{i})}\,,$
(33)

with $\sigma(z)$ representing the error for each SnIa measurement at redshift
$z_{i}$ and where

$\Delta(z_{i})=m^{\text{(obs)}}(z_{i})-5\log{\left(\frac{H_{0}}{c}d_{L}(z_{i})\right)}\,,$
(34)

is the difference between the observed magnitude $m^{\text{(obs)}}$ and the
$H_{0}$-independent luminosity distance.

###  4.2 Cosmic Microwave Background Radiation

In order to assess the behavior of this model at high redshifts, we condense
the information provided by the CMB into what are known as the shift
parameters [33]

$\displaystyle R$ $\displaystyle=$
$\displaystyle\sqrt{\Omega_{b}+\Omega_{c}}\,\frac{H_{0}}{c}\,r(z_{*})\,,$ (35)
$\displaystyle l_{a}$ $\displaystyle=$
$\displaystyle\pi\frac{r(z_{*})}{r_{s}(z_{*})}\,,$ (36)
$\displaystyle\omega_{b}$ $\displaystyle=$ $\displaystyle\Omega_{b}h^{2}\,,$
(37)

where $\Omega_{m}=\Omega_{b}+\Omega_{c}$ is the sum of the baryonic and dark
matter densities, $r(z)$ is the comoving distance

$r(z)=\int_{0}^{z}\frac{c}{H(z)}dz\,,$ (38)

and $r_{s}(z)$ is the comoving sound horizon

$r_{s}(z)=\int_{z}^{\infty}\frac{c_{s}(z)}{H(z)}dz\,,$ (39)

with $c_{s}(z)$, the sound speed as a function of redshift

$c_{s}(z)=\frac{c}{\sqrt{3(1+\hat{R}_{b}/(1+z))}}\,,$ (40)

where $\hat{R}_{b}$ is given by

$\hat{R}_{b}=31500\,\Omega_{b}h^{2}\left(\frac{T_{\rm
CMB}}{2.7K}\right)^{-4}\,,$ (41)

with the temperature of the CMB being taken as $T_{\text{CMB}}=2.7255\,$K.

We can compute the value of the redshift of photon decoupling surface,
$z_{*}$, as [34]

$z_{*}=1048\left[1+\frac{0.00124}{(\Omega_{b}h^{2})^{0.738}}\right]\left[1+g_{1}(\Omega_{m}h^{2})^{g_{2}}\right]\,,$
(42)

where $g_{1}$ and $g_{2}$ are

$\displaystyle g_{1}$ $\displaystyle=$
$\displaystyle\frac{0.0783\,(\Omega_{b}h^{2})^{-0.238}}{1+39.5(\Omega_{b}h^{2})^{-0.763}}\,,$
(43) $\displaystyle g_{2}$ $\displaystyle=$ $\displaystyle
0.560\left[1+21.1\,(\Omega_{b}h^{2})^{1.81}\right]^{-1}\,.$ (44)

Besides the shift parameters, we fix the value of $\omega_{r}=\Omega_{r}h^{2}$
to $\omega_{r}=4.15\times 10^{-5}$.

We consider a Gaussian likelihood of the form

$L=\exp{\left(-\frac{1}{2}\vec{x}\,^{T}C^{-1}\vec{x}\right)}\,,$ (45)

where
$\vec{x}=(l_{a}-l_{a}^{\text{(exp)}},R-R^{\text{(exp)}},\omega_{b}-\omega_{b}^{\text{(exp)}})$
is the vector formed by the theoretical prediction and the measurements of the
three quantities defined above. $C^{-1}$ is the inverse of the covariance
matrix. Both the observational values and the covariance matrix were taken
from the Planck 2018 TT+lowE release [35, 36].

###  4.3 Standard Sirens

Even though this model of $f(Q)$ quickly recovers GR as the redshift
increases, the GW luminosity distance, presented previously in eq. 19, is not
identical to the standard luminosity distance function. Due to the lack of
constraining power coming from current SS events, we instead rely on forecast
data in order to assess the potential that current and future GW observatories
will have in constraining this model.

To generate the mock catalogs of SS events we took the procedure presented in
[37], where mock catalogs for Laser Interferometer Gravitational-Wave
Observatory (LIGO), LISA and the ET were developed. In short, the authors took
the theoretical expected distribution of observed SS events and the expected
error for each measurements for a given observatory, and, assuming a
$\Lambda$CDM cosmological model in agreement with SnIa, generated catalogs of
SS events consisting of a redshift $z$, the value of the luminosity distance
$d_{\text{GW}}(z)$ and its corresponding error $\sigma(z)$.

For this analysis we considered two distinct cases, an ET catalog, with 1000
events, and a LISA catalog, with 15 events and a constraining power
corresponding to the median of the generated catalogs in the previously
mentioned article. We considered mock catalogs for LIGO but did not use them
here as the events provided do not place substantial constraints.

We consider a standard Gaussian likelihood for this dataset of the form

$L=\exp\left(-\frac{1}{2}\sum_{i=0}^{N}\left[\frac{d_{\text{GW}}^{\text{(obs)}}(z_{i})-d_{\text{GW}}(z_{i})}{\sigma(z_{i})}\right]^{2}\right)\,,$
(46)

where $\sigma(z_{i})$ represents the error of the measurement of the SS events
with redshift $z_{i}$ with the corresponding observed luminosity distance
$d_{\text{GW}}^{\text{(obs)}}(z_{i})$ and $N$ represents the number of SS
events in each catalog.

##  5 Constraints

###  5.1 Bayesian Inference Methodology

We performed a Bayesian analysis using Markov chain Monte Carlo (MCMC)
methods. We use PyStan [38], a Python interface to Stan [39], a statistical
programming language that implements the No-U-turn sampler, a variant of the
Hamiltonian Monte Carlo. The output was then analyzed using GetDist [40] and
ArviZ [41] Python packages, that allow for an exploratory analysis of Bayesian
models.

For each run we executed at least four independent chains, each of which with
at least 2500 samples on the posterior distribution and at least 2500 warm-up
steps. To ensure that the chains converged, we used the $\hat{R}$ diagnostic,
a convergence test introduced in [42]. This diagnostic indicates how well the
different chains are mixed together and whether they have converged. We have
only used samples that have $\hat{R}\leq 1.05$, where $\hat{R}=1$ represents
perfect convergence.

The initial values are randomly sampled from a Gaussian distribution with mean
around the expected value and a standard deviation of order 10% of the
corresponding mean for each parameter.

We considered weakly informative priors on each parameter, specifically a
Gaussian distribution centered on the expected value with a standard deviation
two orders of magnitude larger. This ensures a quasi flat prior around the
region of interest.

###  5.2 Results and Discussion

Using the methodology described in the previous subsection, we used the
datasets presented in section 4 to constrain both the $f(Q)$ model under study
and $\Lambda$CDM.

Starting with the constraints on $\Lambda$CDM, the corner plot is presented in
fig. 2, and the best fit values with the corresponding 1$\sigma$ region are
presented in table 2.

Figure 2:  Constraints set by SnIa, CMB shift parameters and forecast SS events for the ET and a median LISA catalog for $\Lambda$CDM.  | SnIa | CMB | LISA | ET
---|---|---|---|---
$\Omega_{m}$ | $0.285\pm 0.012$ | $0.300\pm 0.009$ | $0.304^{+0.035}_{-0.048}$ | $0.290\pm 0.017$
$h$ |  | $0.684\pm 0.007$ | $0.694\pm 0.016$ | $0.699\pm 0.006$
Table 2: Best fit and 1$\sigma$ error set by SnIa, CMB and forecast SS events
for LISA and the ET, for $\Lambda$CDM.

From fig. 2, we can see that in $\Lambda$CDM there is an agreement in the
value of $\Omega_{m}$ for all datasets, with LISA showing higher error bars
when compared to the other datasets. As for the value of $h$ there is a slight
disagreement between the CMB and the SS mock catalogs, given that the mock
catalogs used, constructed in [37], had $h=0.7$ as fiducial value.

Figure 3: Constraints set by SnIa, CMB shift parameters and forecast SS events for the ET and the median case for LISA, for the model of $f(Q)$ presented in eq. 15. | SnIa | CMB | LISA | ET
---|---|---|---|---
$\Omega_{m}$ | $0.335\pm 0.012$ | $0.262\pm 0.008$ | $0.270^{+0.029}_{-0.039}$ | $0.257\pm 0.016$
$h$ |  | $0.732\pm 0.008$ | $0.689\pm 0.009$ | $0.699\pm 0.006$
Table 3: Best fit and 1$\sigma$ error set by SnIa, CMB and forecast SS events
for LISA and the ET, for the model of $f(Q)$ under study.

For the $f(Q)$ model the corner plot is presented in fig. 3 and the best fit
values, with the corresponding 1$\sigma$ region, are presented in table 3. We
can immediately see that there are tensions in both $\Omega_{m}$ and $h$ in
this model.

The best fit value we obtained for $\Omega_{m}$ using SnIa is consistent with
the analysis made in [27] that also includes baryonic acoustic oscillations
and cosmic chronometers. The value $\Omega_{m}\simeq 0.335$ differs from the
one obtained for $\Lambda$CDM, $\Omega_{m}\simeq 0.285$, not surprisingly
given that the low redshift eq. 18 includes a multiplicative factor
$e^{-\lambda}$.

On the other hand, the SnIa value for $\Omega_{m}$ is in tension with the
value we obtain using CMB. This happens because the CMB takes into account a
cumulative effect from high to low redshifts, thus considering both limits eq.
17 and eq. 18.

We observe that the $\Omega_{m}$ values of SnIa and SS are also in tension. It
is true that both observables probe the same low redshift regime, however
their luminosity distances are different due to the modified gravity factor
seen in eq. 19. This extra factor fully explains the tension between the two
measurements, and can be seen as a direct test of the existence of a modified
gravity term. We have verified that the tension in $\Omega_{m}$ between the SS
events and SnIa disappears removing this correction from the luminosity
distance of GW. This analysis shows that comparing future SS events with
standard candles can be decisive in testing modified gravity models using low
redshift observables alone.

As in the case of the SnIa and CMB, we also expect to see a tension arising
between the GW observables and the CMB. This tension, however, depends on the
specific choices for the fiducial values used in the SS mock data. In this
instance, we can see a clear tension between the measurements of $h$, but
comparable values for $\Omega_{m}$. This latter result is just a coincidence
occurring for this particular choice of $h$, and is unlikely to be a general
feature.

##  6 Conclusions

We studied an $f(Q)$ modification of gravity that accounts for the late time
accelerated expansion of the Universe without adding a cosmological constant.
This modification has the same number of free parameters as the standard
$\Lambda$CDM model, albeit, the propagation of tensorial perturbations
introduces a change in the luminosity distance of GWs.

We used a dynamical system analysis to reveal the existence of three fixed
points and three disjoint regions in the phase space defined by the sign of
the parameter $\lambda$. For $\lambda<0$ we identified a region characterized
by trajectories where $E(z)\rightarrow 0$ as the scale factor increases. When
$\lambda=0$, we recover a CDM universe with a radiation followed by a matter
epoch. When $\lambda>0$, the stable fixed point corresponds to a regime of
eternal inflation, meaning that this case can lead to a viable universe with a
late-time accelerated expansion without an explicit cosmological constant.

First we constrained this model using SnIa recovering the value of
$\Omega_{m}$ obtained previously in [27], where other low redshift
observations were also used. This value is, not surprisingly, different from
the best fit found for a $\Lambda$CDM, given the factor
$e^{-\lambda}\Omega_{m}$ in the modified Friedmann equation eq. 18.

We made an independent constraint using high redshift measurements from CMB
and verified that the best fit value for $\Omega_{m}$ is in tension with the
value obtained with SnIa. We conclude that this is a manifestation of
differing high and low redshift limits of the model, that CMB takes into
account but SnIa does not. In order to keep the attractive features of this
model at low redshifts, the form of $f(Q)$ needs to be somehow modified at
high redshifts to make it consistent with CMB data.

As is, the current data disfavors this model already at the background level.
This is a common feature in models of modified gravity that depart from
$\Lambda$CDM. $f(T)$ and $f(Q)$ models can be shown to be equivalent in FLRW,
whereas $f(R)$ models, though being formally distinguishable, must in practice
have a similar cosmological evolution. In particular, $f(R)$ models have their
free parameters strongly constrained to be very close to the $\Lambda$CDM
limit [43], which seems also to be the case for $f(T)$ and $f(Q)$ [44, 45,
16].

Finally we used mock catalogs of SS constructed with $\Lambda$CDM as the
fiducial model. As expected, the modified luminosity distance for
gravitational waves raises a tension with the SnIa results implying that in
the future, this feature will immediately discriminate this model from
$\Lambda$CDM.

Although current data (CMB and SnIa) already disfavor this model, this study
provides a working example of how future SS events, measured both by LISA and
the ET, will play a fundamental role when testing models of modified gravity
using low redshift data alone.

###### Acknowledgements.

The authors are thankful to Fotis Anagnostopoulos for the fruitful
discussions. The authors acknowledge support from Fundação para a Ciência e a
Tecnologia via the following projects: UIDB/04434/2020 & UIDP/04434/2020,
CERN/FIS-PAR/0037/2019, PTDC/FIS-OUT/29048/2017, COMPETE2020:
POCI-01-0145-FEDER-028987 & FCT: PTDC/FIS-AST/28987/2017, PTDC/FIS-
AST/0054/2021, EXPL/FIS-AST/1368/2021.

## References

  * Perivolaropoulos and Skara [2022] L. Perivolaropoulos and F. Skara, Challenges for ${\Lambda}$CDM: An update, New Astronomy Reviews 95, 101659 (2022), arXiv:2105.05208 [astro-ph.CO] .
  * Jiménez  _et al._ [2019] J. B. Jiménez, L. Heisenberg, and T. Koivisto, The geometrical trinity of gravity, Universe 5, 173 (2019), arXiv:1903.06830 [hep-th] .
  * Jimenez  _et al._ [2017] J. B. Jimenez, L. Heisenberg, and T. Koivisto, Coincident General Relativity, Phys. Rev. D 98, 044048 (2018) 10.1103/PhysRevD.98.044048 (2017), arXiv:1710.03116 [gr-qc] .
  * Cai and Yang [2016] R.-G. Cai and T. Yang, Estimating cosmological parameters by the simulated data of gravitational waves from the Einstein Telescope, Phys. Rev. D 95, 044024 (2017) 10.1103/PhysRevD.95.044024 (2016), arXiv:1608.08008 [astro-ph.CO] .
  * Saridakis  _et al._ [2021] E. N. Saridakis, R. Lazkoz, V. Salzano, P. V. Moniz, S. Capozziello, J. B. Jiménez, M. D. Laurentis, and G. J. Olmo, eds.,  _Modified Gravity and Cosmology_ (Springer International Publishing, 2021).
  * Ayuso  _et al._ [2022] I. Ayuso, R. Lazkoz, and J. P. Mimoso, DGP and DGP-like cosmologies from $f({Q})$ actions, Phys. Rev. D 105, 083534 (2022).
  * Barros  _et al._ [2020] B. J. Barros, T. Barreiro, T. Koivisto, and N. J. Nunes, Testing $F({Q})$ gravity with redshift space distortions, Physics of the Dark Universe 30, 100616 (2020), arXiv:2004.07867 [gr-qc] .
  * Belgacem  _et al._ [2017] E. Belgacem, Y. Dirian, S. Foffa, and M. Maggiore, The gravitational-wave luminosity distance in modified gravity theories, Phys. Rev. D 97, 104066 (2018) 10.1103/PhysRevD.97.104066 (2017), arXiv:1712.08108 [astro-ph.CO] .
  * Albuquerque and Frusciante [2022] I. S. Albuquerque and N. Frusciante, A designer approach to f(Q) gravity and cosmological implications, Physics of the Dark Universe 35, 100980 (2022).
  * Atayde and Frusciante [2023] L. Atayde and N. Frusciante, f(Q) gravity and neutrino physics, Phys. Rev. D 107, 124048 (2023), arXiv:2306.03015 [astro-ph.CO] .
  * Atayde and Frusciante [2021] L. Atayde and N. Frusciante, Can $f({Q})$ gravity challenge $\mathrm{\Lambda}\mathrm{CDM}$?, Phys. Rev. D 104, 064052 (2021).
  * Frusciante [2021] N. Frusciante, Signatures of $f({Q})$ gravity in cosmology, Phys. Rev. D 103, 044021 (2021).
  * Capozziello  _et al._ [2022] S. Capozziello, V. D. Falco, and C. Ferrara, Comparing equivalent gravities: common features and differences, The European Physical Journal C 82, 10.1140/epjc/s10052-022-10823-x (2022).
  * Krishnan  _et al._ [2021] C. Krishnan, E. Ó. Colgáin, M. Sheikh-Jabbari, and T. Yang, Running hubble tension and a h0 diagnostic, Physical Review D 103, 10.1103/physrevd.103.103509 (2021).
  * Khyllep  _et al._ [2021] W. Khyllep, A. Paliathanasis, and J. Dutta, Cosmological solutions and growth index of matter perturbations in $f(Q)$ gravity, Phys. Rev. D 103, 103521 (2021), arXiv:2103.08372 [gr-qc] .
  * Ayuso  _et al._ [2021] I. Ayuso, R. Lazkoz, and V. Salzano, Observational constraints on cosmological solutions of $f({Q})$ theories, Phys. Rev. D 103, 063505 (2021), 2012.00046 [astro-ph.CO] .
  * Collaboration and Collaboration [2017] T. L. S. Collaboration and T. V. Collaboration, Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119 161101 (2017) 10.1103/PhysRevLett.119.161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Collaboration and the Virgo Collaboration [2020] T. L. S. Collaboration and the Virgo Collaboration, Properties and astrophysical implications of the 150 msun binary black hole merger GW190521, Astrophys. J. Lett. 900, L13 (2020) 10.3847/2041-8213/aba493 (2020), arXiv:2009.01190 [astro-ph.HE] .
  * Graham  _et al._ [2020] M. J. Graham, K. E. S. Ford, B. McKernan, N. P. Ross, D. Stern, K. Burdge, M. Coughlin, S. G. Djorgovski, A. J. Drake, D. Duev, M. Kasliwal, A. A. Mahabal, S. van Velzen, J. Belicki, E. C. Bellm, R. Burruss, S. B. Cenko, V. Cunningham, G. Helou, S. R. Kulkarni, F. J. Masci, T. Prince, D. Reiley, H. Rodriguez, B. Rusholme, R. M. Smith, and M. T. Soumagnac, Candidate electromagnetic counterpart to the binary black hole merger gravitational wave event S190521g, Phys. Rev. Lett. 124, 251102 10.1103/PhysRevLett.124.251102 (2020), arXiv:2006.14122 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, J. Camp, C. Caprini, V. Cardoso, M. Colpi, J. Conklin, N. Cornish, C. Cutler, K. Danzmann, R. Dolesi, L. Ferraioli, V. Ferroni, E. Fitzsimons, J. Gair, L. G. Bote, D. Giardini, F. Gibert, C. Grimani, H. Halloin, G. Heinzel, T. Hertog, M. Hewitson, K. Holley-Bockelmann, D. Hollington, M. Hueller, H. Inchauspe, P. Jetzer, N. Karnesis, C. Killow, A. Klein, B. Klipstein, N. Korsakova, S. L. Larson, J. Livas, I. Lloro, N. Man, D. Mance, J. Martino, I. Mateos, K. McKenzie, S. T. McWilliams, C. Miller, G. Mueller, G. Nardini, G. Nelemans, M. Nofrarias, A. Petiteau, P. Pivato, E. Plagnol, E. Porter, J. Reiche, D. Robertson, N. Robertson, E. Rossi, G. Russano, B. Schutz, A. Sesana, D. Shoemaker, J. Slutsky, C. F. Sopuerta, T. Sumner, N. Tamanini, I. Thorpe, M. Troebs, M. Vallisneri, A. Vecchio, D. Vetrugno, S. Vitale, M. Volonteri, G. Wanner, H. Ward, P. Wass, W. Weber, J. Ziemer, and P. Zweifel, Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Auclair  _et al._ [2023] P. Auclair, D. Bacon, T. Baker, T. Barreiro, N. Bartolo, E. Belgacem, N. Bellomo, I. Ben-Dayan, D. Bertacca, M. Besancon, J. J. Blanco-Pillado, D. Blas, G. Boileau, G. Calcagni, R. Caldwell, C. Caprini, C. Carbone, C.-F. Chang, H.-Y. Chen, N. Christensen, S. Clesse, D. Comelli, G. Congedo, C. Contaldi, M. Crisostomi, D. Croon, Y. Cui, G. Cusin, D. Cutting, C. Dalang, V. De Luca, W. D. Pozzo, V. Desjacques, E. Dimastrogiovanni, G. C. Dorsch, J. M. Ezquiaga, M. Fasiello, D. G. Figueroa, R. Flauger, G. Franciolini, N. Frusciante, J. Fumagalli, J. García-Bellido, O. Gould, D. Holz, L. Iacconi, R. K. Jain, A. C. Jenkins, R. Jinno, C. Joana, N. Karnesis, T. Konstandin, K. Koyama, J. Kozaczuk, S. Kuroyanagi, D. Laghi, M. Lewicki, L. Lombriser, E. Madge, M. Maggiore, A. Malhotra, M. Mancarella, V. Mandic, A. Mangiagli, S. Matarrese, A. Mazumdar, S. Mukherjee, I. Musco, G. Nardini, J. M. No, T. Papanikolaou, M. Peloso, M. Pieroni, L. Pilo, A. Raccanelli, S. Renaux-Petel, A. I. Renzini, A. Ricciardone, A. Riotto, J. D. Romano, R. Rollo, A. R. Pol, E. R. Morales, M. Sakellariadou, I. D. Saltas, M. Scalisi, K. Schmitz, P. Schwaller, O. Sergijenko, G. Servant, P. Simakachorn, L. Sorbo, L. Sousa, L. Speri, D. A. Steer, N. Tamanini, G. Tasinato, J. Torrado, C. Unal, V. Vennin, D. Vernieri, F. Vernizzi, M. Volonteri, J. M. Wachter, D. Wands, L. T. Witkowski, M. Zumalacárregui, J. Annis, F. R. Ares, P. P. Avelino, A. Avgoustidis, E. Barausse, A. Bonilla, C. Bonvin, P. Bosso, M. Calabrese, M. Çalışkan, J. A. R. Cembranos, M. Chala, D. Chernoff, K. Clough, A. Criswell, S. Das, A. d. Silva, P. Dayal, V. Domcke, R. Durrer, R. Easther, S. Escoffier, S. Ferrans, C. Fryer, J. Gair, C. Gordon, M. Hendry, M. Hindmarsh, D. C. Hooper, E. Kajfasz, J. Kopp, S. M. Koushiappas, U. Kumar, M. Kunz, M. Lagos, M. Lilley, J. Lizarraga, F. S. N. Lobo, A. Maleknejad, C. J. A. P. Martins, P. D. Meerburg, R. Meyer, J. P. Mimoso, S. Nesseris, N. Nunes, V. Oikonomou, G. Orlando, O. Özsoy, F. Pacucci, A. Palmese, A. Petiteau, L. Pinol, S. P. Zwart, G. Pratten, T. Prokopec, J. Quenby, S. Rastgoo, D. Roest, K. Rummukainen, C. Schimd, A. Secroun, A. Sesana, C. F. Sopuerta, I. Tereno, A. Tolley, J. Urrestilla, E. C. Vagenas, J. van de Vis, R. van de Weygaert, B. Wardell, D. J. Weir, G. White, B. Świeżewska, V. I. Zhdanov, and T. L. C. W. Group, Cosmology with the laser interferometer space antenna, Living Reviews in Relativity 26, 5 (2023).
  * [22] https://github.com/jpmvferreira/testing-L-free-fQ-cosmology.
  * Järv  _et al._ [2018] L. Järv, M. Rünkla, M. Saal, and O. Vilson, Nonmetricity formulation of general relativity and its scalar-tensor extension, Phys. Rev. D 97, 124025 (2018), arXiv:1802.00492 [gr-qc] .
  * Jiménez  _et al._ [2019] J. B. Jiménez, L. Heisenberg, T. S. Koivisto, and S. Pekar, Cosmology in $f({Q})$ geometry, Phys. Rev. D 101, 103507 (2020) 10.1103/PhysRevD.101.103507 (2019), arXiv:1906.10027 [gr-qc] .
  * Jiménez  _et al._ [2018] J. Jiménez, L. Heisenberg, and T. S. Koivisto, Teleparallel palatini theories, Journal of Cosmology and Astroparticle Physics 08 (08), 039, arXiv:1803.10185 [gr-qc] .
  * Koivisto [2006] T. Koivisto, A note on covariant conservation of energy–momentum in modified gravities, Classical and Quantum Gravity 23, 4289 (2006).
  * Anagnostopoulos  _et al._ [2021] F. K. Anagnostopoulos, S. Basilakos, and E. N. Saridakis, First evidence that non-metricity f(Q) gravity could challenge ${\Lambda}$CDM, Physics Letters B 822, 136634 (2021).
  * Anagnostopoulos  _et al._ [2023] F. K. Anagnostopoulos, V. Gakis, E. N. Saridakis, and S. Basilakos, New models and big bang nucleosynthesis constraints in f(Q) gravity, The European Physical Journal C 83, 58 (2023).
  * Bahamonde  _et al._ [2018] S. Bahamonde, C. G. Böhmer, S. Carloni, E. J. Copeland, W. Fang, and N. Tamanini, Dynamical systems applied to cosmology: Dark energy and modified gravity, Physics Reports 775-777, 1 (2018), dynamical systems applied to cosmology: Dark energy and modified gravity.
  * Scolnic  _et al._ [2018] D. M. Scolnic, D. O. Jones, A. Rest, Y. C. Pan, R. Chornock, R. J. Foley, M. E. Huber, R. Kessler, G. Narayan, A. G. Riess, S. Rodney, E. Berger, D. J. Brout, P. J. Challis, M. Drout, D. Finkbeiner, R. Lunnan, R. P. Kirshner, N. E. Sanders, E. Schlafly, S. Smartt, C. W. Stubbs, J. Tonry, W. M. Wood-Vasey, M. Foley, J. Hand, E. Johnson, W. S. Burgett, K. C. Chambers, P. W. Draper, K. W. Hodapp, N. Kaiser, R. P. Kudritzki, E. A. Magnier, N. Metcalfe, F. Bresolin, E. Gall, R. Kotak, M. McCrum, and K. W. Smith (Pan-STARRS1), The complete light-curve sample of spectroscopically confirmed Type Ia Supernovae from Pan-STARRS1 and cosmological constraints from the combined Pantheon sample, Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .
  * [31] https://github.com/dscolnic/Pantheon.
  * Goliath  _et al._ [2001] M. Goliath, R. Amanullah, P. Astier, A. Goobar, and R. Pain, Supernovae and the nature of the dark energy, Astronomy & Astrophysics 380, 6 (2001), arXiv:astro-ph/0104009 [astro-ph] .
  * Zhai and Wang [2019] Z. Zhai and Y. Wang, Robust and model-independent cosmological constraints from distance measurements, JCAP 07, 005, arXiv:1811.07425 [astro-ph.CO] .
  * Hu and Sugiyama [1996] W. Hu and N. Sugiyama, Small-scale cosmological perturbations: An analytic approach, The Astrophysical Journal 471, 542 (1996).
  * Perez  _et al._ [2023] J. d. C. Perez, J. S. Peracaula, and C. P. Singh, Running vacuum in brans-dicke theory: a possible cure for the $\sigma_{8}$ and ${H}_{0}$ tensions, (2023), arXiv:2302.04807 [astro-ph.CO] .
  * Collaboration [2020] P. Collaboration, Planck 2018 results, Astronomy & Astrophysics 641, A6 (2020).
  * Ferreira  _et al._ [2022] J. Ferreira, T. Barreiro, J. Mimoso, and N. J. Nunes, Forecasting F(Q) cosmology with $\Lambda$CDM background using standard sirens, Phys. Rev. D 105, 123531 (2022), arXiv:2203.13788 [astro-ph.CO] .
  * Riddell  _et al._ [2021] A. Riddell, A. Hartikainen, and M. Carter, pystan (3.0.0), PyPi (2021).
  * Team [2021] S. D. Team,  _Stan Modeling Language Users Guide and Reference Manual_, 2nd ed. (2021).
  * Lewis [2019] A. Lewis, Getdist: a python package for analysing monte carlo samples, (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Kumar  _et al._ [2019] R. Kumar, C. Carroll, A. Hartikainen, and O. Martin, ArviZ a unified library for exploratory analysis of bayesian models in python, Journal of Open Source Software 4, 1143 (2019).
  * Vehtari  _et al._ [2021] A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, and P.-C. Bürkner, Rank-normalization, folding, and localization: An improved Rˆ for assessing convergence of MCMC (with discussion), Bayesian Analysis 16, 10.1214/20-ba1221 (2021).
  * Nunes  _et al._ [2017] R. C. Nunes, S. Pan, E. N. Saridakis, and E. M. Abreu, New observational constraints on $f({R})$ gravity from cosmic chronometers, Journal of Cosmology and Astroparticle Physics 2017 (01), 005.
  * dos Santos  _et al._ [2022] F. B. M. dos Santos, J. E. Gonzalez, and R. Silva, Observational constraints on f(T) gravity from model-independent data, The European Physical Journal C 82, 10.1140/epjc/s10052-022-10784-1 (2022).
  * Lazkoz  _et al._ [2019] R. Lazkoz, F. S. N. Lobo, M. Ortiz-Baños, and V. Salzano, Observational constraints of $f({Q})$ gravity, Physical Review D 100, 104027 (2019), arXiv:1907.13219 [gr-qc] .