###### Abstract

The early Cosmology driven by a single scalar field, both massless and
massive, in the context of Eddington-inspired Born-Infeld gravity, is
explored. We show the existence of nonsingular solutions of bouncing and
loitering type (depending on the sign of the gravitational theory’s parameter,
$\epsilon$) replacing the Big Bang singularity, and discuss their properties.
In addition, in the massive case we find some new features of the cosmological
evolution depending on the value of the mass parameter, including asymmetries
in the expansion/contraction phases, or a continuous transition between a
contracting phase to an expanding one via an intermediate loitering phase. We
also provide a combined analysis of cosmic chronometers, standard candles,
BAO, and CMB data to constrain the model, finding that for roughly
$|\epsilon|\lesssim 5\cdot 10^{-8}\text{ m}^{2}$ the model is compatible with
the latest observations while successfully removing the Big Bang singularity.
This bound is several orders of magnitude stronger than the most stringent
constraints currently available in the literature.

xx 1 5 Received: date; Accepted: date; Published: date Singularity-free and
cosmologically viable Born-Infeld gravity with scalar matter David Benisty1,2,
Gonzalo J. Olmo3, Diego Rubiera-Garcia∗4 Correspondence<EMAIL_ADDRESS>
##  1 Introduction

The standard concordance cosmological $\Lambda$CDM model, framed within
Einstein’s General Theory of Relativity (GR), including an early phase of
inflationary expansion, a cold dark matter component, and a tiny cosmological
constant driving the accelerated late-time expansion of the Universe, has
successfully met all observations Amendola et al. (2013); Aghanim et al.
(2020). Within this model, scalar fields have found new and imaginative
applications. For instance, inflationary models in the early Universe involve
from one to many scalar fields Starobinsky (1979, 1980); Guth (1981); Albrecht
and Steinhardt (1982); Mukhanov and Chibisov (1981); Guth and Pi (1982); Linde
(1982); Barrow and Cotsakis (1988); Barrow (1988); Elizalde et al. (2008). In
the slow-roll approximation the exact form of the scalar field potential is
unknown since many different potentials have been studied and confronted to
observations. On the other hand, a way to parameterize dark energy is by using
a scalar field, the so-called quintessence model (for canonical scalar fields
Ratra and Peebles (1988); Caldwell et al. (1998)) or its generalizations to
K-essence models (when a non-canonical scalar Lagrangian is considered
Kehayias and Scherrer (2019); Oikonomou and Chatzarakis (2020); Chakraborty et
al. (2019); Babichev et al. (2018)), in such a way that the cosmological
constant gets replaced by a dark energy fluid with a nearly constant density
today Zlatev et al. (1999); Caldwell (2002); Chiba et al. (2000); Bento et al.
(2002); Tsujikawa (2013). Dark matter can be also parameterized in terms of
weakly-interacting massive particles, which can be scalar particles still
undiscovered at colliders and other dark matter detection experiments. Models
for dark matter can also be based on other kinds of scalar fields, for
instance, via fuzzy dark matter Hu et al. (2000), or by using a Lagrange
multiplier that changes the behaviour of the kinetic term Anagnostopoulos et
al. (2019); Benisty et al. (2019); Benisty and Guendelman (2018, 2017). Scalar
field models may also be the result of complex effective interactions of other
fundamental fields in equilibrium, such as in Bose-Einstein condensates, thus
allowing for an even broader range of phenomenological justifications.
Therefore, the consideration of scalar fields within cosmological models is
justified from the point of view of building a fully consistent history of the
cosmological evolution compatible with observations.

On the gravity side, and despite its observational success, from a fundamental
point of view the $\Lambda$CDM model still contains a visible singularity at
the Universe’s past. This is an unavoidable consequence of the singularity
theorems (see e.g. Senovilla and Garfinkle (2015) for a pedagogical
discussion). In order to tackle this issue, it is widely assumed in the
community that at the strong curvatures and fields of the very early Universe,
quantum gravity effects should come into play in order to regularize this
singularity. Since a quantum theory of gravity is not available yet, an
effective way to capture such hypothetical effects is via modified theories of
gravity De Felice and Tsujikawa (2010); Capozziello and De Laurentis (2011);
Nojiri et al. (2017); Heisenberg (2019). The plausibility of such alternative
gravitational descriptions to supersede GR and represent observationally
viable alternatives to the $\Lambda$CDM model has been widely discussed in the
literature according to different perspectives and approaches Bull et al.
(2016).

Among the large pool of theories which have been investigated in the
literature, for the sake of this paper we bring here the proposal originally
introduced by Banados and Ferreira Banados and Ferreira (2010) and dubbed as
Eddington-inspired Born-Infeld (EiBI) gravity111This proposal can actually be
framed within the tradition of considering square-root action, such as in the
DBI one, see e.g. Alishahiha et al. (2004); Liu et al. (2012); Choudhury and
Pal (2013, 2015).. In order to avoid troubles with ghost-like instabilities,
this theory is typically formulated in metric-affine spaces, where metric and
affine connection are a priori independent entities. EiBI theory has found
many different applications in astrophysics and cosmology, see e.g. Harko et
al. (2013); Wei et al. (2015); Shaikh (2015); Avelino (2016); Prasetyo et al.
(2018); Chen et al. (2018); Shaikh (2018); Jana et al. (2018); Böhmer and
Fiorini (2019); Delhom et al. (2019). In particular, the existence of bouncing
solutions replacing the Big Bang singularity within EiBI gravity was first
hinted in Avelino and Ferreira (2012), in particular, when scalar fields are
considered as the matter source. Further works in the subject in the last few
years have reinforced the capability of this theory to remove such
singularities according to different mechanisms (for a review, see Beltran
Jimenez et al. (2018)).

The main aim of this work is to construct explicit such singularity-free
solutions within the early cosmological evolution, corresponding to a single
massless (quintaessential) scalar field, and to further extend (numerically)
this analysis to the massive case. We shall show the existence of two kinds of
singularity-free solutions depending on the sign of the EiBI gravity
parameter. The first one corresponds to bouncing solutions, where the universe
contracts down to a minimum size before entering into an expansion phase,
while the second are loitering solutions, which interpolate between an
asymptotically Minkowski past and the current cosmological evolution. For
these nonsingular solutions we carry out a combined analysis of cosmic
chronometers, standard candles, BAO, and CMB data in order to constrain the
EiBI parameter, finding the bound $|\epsilon|\lesssim 5\cdot 10^{-8}\text{
m}^{2}$.

This paper is organized as follows: in Sec.2 we introduce EiBI gravity,
discuss its properties, and construct its cosmological equations when coupled
to a scalar field. In Sec.3 we consider massless scalar fields and discuss its
corresponding bouncing and loitering solutions, extending these results in
Sec.4 to the massive case. Sec.5 sees the theory fit with the latest
observations (on $\Lambda$CDM background), and we conclude in Sec.6 with a
summary and some perspectives.

##  2 Eddington-inspired Born-Infeld gravity

###  2.1 Action and basic field equations

The action of EiBI gravity can be conveniently written as

$\mathcal{S}_{EiBI}=\frac{1}{\epsilon\kappa^{2}}\int
d^{4}x\left[\sqrt{-q}-\lambda\sqrt{-g}\right]+\mathcal{S}_{m}(g_{\mu\nu},\psi_{m})\
,$ (1)

where $\kappa^{2}\equiv 8\pi G/c^{4}$ is Newton’s constant, $\epsilon$ is EiBI
parameter, $\lambda$ is a dimensionless constant, $g$ is the determinant of
the space-time metric $g_{\mu\nu}$ and $q$ the determinant of an auxiliary
metric defined as:

$q_{\mu\nu}\equiv g_{\mu\nu}+\epsilon R_{(\mu\nu)}(\Gamma)\ ,$ (2)

where the (symmetric part of the) Ricci tensor
$R_{\mu\nu}(\Gamma)\equiv{R^{\rho}}_{\mu\rho\nu}(\Gamma)$ is a function solely
of the (torsionless) affine connection
$\Gamma\equiv\Gamma_{\mu\nu}^{\lambda}$, assumed to be a priori independent of
the space-time metric $g_{\mu\nu}$ (metric-affine or Palatini formalism). This
symmetrizing requirement ensures that the theory is invariant under projective
transformations, which avoids the presence of ghost-like instabilities Beltrán
Jiménez and Delhom (2019, 2020). As for the matter sector,
$\mathcal{S}_{m}=\int d^{4}x\sqrt{-g}\mathcal{L}_{m}(g_{\mu\nu},\psi_{m})$, it
is assumed to be minimally coupled to the space-time metric $g_{\mu\nu}$ (see
e.g. Delhom (2020) for a definition of minimal coupling in metric-affine
theories), with $\psi_{m}$ denoting collectively the matter fields. It is
worth pointing out that EiBI gravity recovers the GR dynamics and its
solutions in the $|R_{\mu\nu}|\ll\epsilon^{-1}$ limit Beltran Jimenez et al.
(2018) (from now on, vertical bars indicate a determinant), inheriting an
effective cosmological constant
$\Lambda_{eff}=\frac{\lambda-1}{\epsilon\kappa^{2}}$. This fact allows EiBI
gravity to naturally pass weak-field limit tests, such as those based on solar
system observations.

EiBI gravity is nowadays a very well known theory in the community thanks to
its many applications (for a detailed account of its properties we refer the
reader to the review Beltran Jimenez et al. (2018)). It is actually a member
of the so-called Ricci-based family of gravitational theories Afonso et al.
(2018), all of which admit an Einstein-like representation of their field
equations given by

${G^{\mu}}_{\nu}(q)=\frac{\kappa^{2}}{|\Omega|^{1/2}}\left[{T^{\mu}}_{\nu}-\delta^{\mu}_{\nu}\left(\mathcal{L}_{G}+\frac{T}{2}\right)\right]\
,$ (3)

where $\mathcal{L}_{G}$ is the gravitational Lagrangian,
$T_{\mu\nu}=\frac{2}{\sqrt{-g}}\frac{\delta S_{m}}{\delta g^{\mu\nu}}$ is the
stress-energy tensor of the matter fields, $T$ its trace, and
${G^{\mu}}_{\nu}(q)$ is the Einstein tensor of the auxiliary metric
$q_{\mu\nu}$, such that $\Gamma$ is Levi-Civita of it, that is

$\nabla_{\alpha}(\sqrt{-q}\,q^{\mu\nu})=0\ .$ (4)

This $q_{\mu\nu}$ metric is related to the space-time one $g_{\mu\nu}$ via the
relation

$q_{\mu\nu}=g_{\mu\alpha}\,{\Omega^{\alpha}}_{\nu}\ ,$ (5)

where the _deformation_ matrix $\hat{\Omega}$ depends on-shell on the matter
fields (and possibly the space-time metric $g_{\mu\nu}$ as well). For EiBI
gravity this matrix is implicitly determined via the equation

$|\Omega|^{1/2}{(\Omega^{-1})^{\mu}}_{\nu}=\lambda\delta^{\mu}_{\nu}-\epsilon\kappa^{2}{T^{\mu}}_{\nu}\
,$ (6)

while the EiBI Lagrangian in (1) can also be expressed in terms of this matrix
as

$\mathcal{L}_{G}=\frac{|\Omega|^{1/2}-\lambda}{\epsilon\kappa^{2}}\ .$ (7)

Let us point out that all terms on the right-hand side of the equations (3)
are functions of the matter fields and the metric $g_{\mu\nu}$, thus
representing a system of second-order field equations with new couplings
engendered by the matter fields. In vacuum, ${T^{\mu}}_{\nu}=0$, one recovers
the GR solutions222As a remark, we would like to note that it has been shown
that in generic RBGs, due to the non-linearities of the field equations, the
deformation matrix admits other (typically pathological) solutions besided the
one that boils down to GR in vacuum. These solutions can lead to anisotropic
deformation matrices even if the stress-energy tensor is isotropic.
Remarkably, it was shown that for EiBI gravity, no such anisotropic solutions
exist in presence of an isotropic stress-energy tensor, see Beltrán Jiménez et
al. (2020)., which ensures the propagation of the two polarizations of the
gravitational field travelling at the speed of light Jana et al. (2018).

###  2.2 EiBI cosmology with scalar fields

As the matter sector of our model let us consider a single (real) scalar field
described by the action

$\mathcal{S}_{m}=-\frac{1}{2}\int
d^{4}x\sqrt{-g}\mathcal{L}_{m}=-\frac{1}{2}\int d^{4}x\sqrt{-g}(X+2V(\phi))\
,$ (8)

with kinetic term $X\equiv g^{\mu\nu}\partial_{\mu}\phi\,\partial_{\nu}\phi$
and scalar potential $V(\phi)$. The variation of this action with respect to
the scalar field leads to the field equations

$\frac{1}{\sqrt{-g}}\partial_{\mu}\left(\sqrt{-g}g^{\mu\nu}\partial_{\nu}\phi\right)=V_{\phi}\
,$ (9)

(where $V_{\phi}\equiv dV/d\phi$) and to the associated stress-energy tensor

${T^{\mu}}_{\nu}=g^{\mu\alpha}\partial_{\alpha}\phi\partial_{\nu}\phi-\frac{\mathcal{L}_{m}}{2}\delta^{\mu}_{\nu}\
.$ (10)

We consider next a spatially flat Friedman-Lemaitre-Robertson-Walker (FLRW)
space-time333For simplicity in this work we do not consider the cases of open
and closed universes, though the analysis could also be carried out in such
cases following, for instance, the results of Barragan and Olmo (2010)., with
line element:

$ds_{g}^{2}\equiv g_{\mu\nu}dx^{\mu}dx^{\nu}=-dt^{2}+a^{2}(t)d\vec{x}^{2}\ ,$
(11)

where $a(t)$ is the expansion factor and
$d\vec{x}^{\,2}=\delta_{ij}\,dx^{i}dx^{j}$, with $i,j=1\ldots 3$ for the
spatial part. As from now on we shall assume $\phi=\phi(t)$, the stress-energy
tensor (10) reads:

${T^{\mu}}_{\nu}=\left(\begin{array}[]{cc}-\tfrac{1}{2}\dot{\phi}^{2}-V&0\\\
0&(\tfrac{1}{2}\dot{\phi}^{2}-V){\bf\it I}_{3\times 3}\end{array}\right)\ .$
(12)

where a dot denotes a derivative with respect to $t$, as usual.

From Eq.(6) we can conclude that the deformation matrix needs to have a
similar algebraic (diagonal) structure as that of the stress-energy tensor of
the scalar field, namely

${\Omega^{\mu}}_{\nu}=\left(\begin{array}[]{cc}\Omega_{+}&0\\\
0&\Omega_{-}\,{\bf\it I}_{3\times 3}\end{array}\right)\ ,$ (13)

where the identity matrix $I_{3\times 3}$ represents the spatial sector. Now,
plugging this ansatz and Eq.(12) into Eq.(6) gives the components of this
matrix as:

$\displaystyle\Omega_{-}^{2}$ $\displaystyle=$
$\displaystyle\left(\lambda+\epsilon\kappa^{2}V-\tfrac{\epsilon\kappa^{2}\dot{\phi}^{2}}{2}\right)\left(\lambda+\epsilon\kappa^{2}V+\tfrac{\epsilon\kappa^{2}\dot{\phi}^{2}}{2}\right)$
(14) $\displaystyle=$ $\displaystyle\left(\tilde{\lambda}^{2}-\Phi^{2}\right)$
$\displaystyle\Omega_{+}^{2}$ $\displaystyle=$
$\displaystyle\frac{\left(\lambda+\epsilon\kappa^{2}V-\tfrac{\epsilon\kappa^{2}\dot{\phi}^{2}}{2}\right)^{3}}{\left(\lambda+\epsilon\kappa^{2}V+\tfrac{\epsilon\kappa^{2}\dot{\phi}^{2}}{2}\right)}=\frac{\left(\tilde{\lambda}-\Phi\right)^{3}}{\left(\tilde{\lambda}+\Phi\right)}\
,$ (15)

where we have used the shorthand notations

$\displaystyle\tilde{\lambda}$ $\displaystyle\equiv$
$\displaystyle\lambda+\epsilon\kappa^{2}V$ (16) $\displaystyle\Phi$
$\displaystyle\equiv$
$\displaystyle\tfrac{\epsilon\kappa^{2}\dot{\phi}^{2}}{2}\ .$ (17)

Therefore, we are ready to cast the field equations (3) for this problem as

$\epsilon{R^{\mu}}_{\nu}(q)=\left(\begin{array}[]{cc}1-\frac{\left(\tilde{\lambda}+\Phi\right)}{|\Omega|^{\frac{1}{2}}}&0\\\
0&\left(1-\frac{\left(\tilde{\lambda}-\Phi\right)}{|\Omega|^{\frac{1}{2}}}\right)\delta^{i}_{j}\end{array}\right)\
,$ (18)

where the determinant of the deformation matrix reads
$|\Omega|^{\frac{1}{2}}=\Omega_{+}^{1/2}\Omega_{-}^{3/2}$. Now, to compute the
left-hand side of the field equations (18) one can write another line element
for the auxiliary metric also of the FLRW form, that is

$\displaystyle ds_{q}^{2}$ $\displaystyle\equiv$ $\displaystyle
q_{\mu\nu}dx^{\mu}dx^{\nu}=-dT^{2}+\tilde{a}^{2}(T)d\vec{x}^{2}$
$\displaystyle=$
$\displaystyle-\Omega_{+}(t)dt^{2}+a^{2}(t)\Omega_{-}(t)\,d\vec{x}^{2}\ ,$

where in the second line we have used the fundamental relation (5) together
with the ansatz (13). These gravitational field equations must be supplemented
with the scalar field equations (9) which, in the FLRW background (11), read

$\ddot{\phi}+3\frac{\dot{a}}{a}\dot{\phi}+V_{\phi}=0\ .$ (20)

Now, using the fact that in the $q_{\mu\nu}$ geometry (2.2) we have the well
known formulas

$\displaystyle{R^{t}}_{t}$ $\displaystyle\equiv$
$\displaystyle\frac{3}{\tilde{a}}\frac{d^{2}\tilde{a}}{dT^{2}}$ (21)
$\displaystyle{R^{i}}_{i}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\tilde{a}}\frac{d^{2}\tilde{a}}{dT^{2}}+\frac{2}{\tilde{a}^{2}}\left(\frac{d\tilde{a}}{dT}\right)^{2}\,,$
(22)

with the relations

$\displaystyle\frac{d\tilde{a}}{dT}$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega_{+}^{\frac{1}{2}}}\frac{d}{dt}\left(a\Omega_{-}^{\frac{1}{2}}\right)$
(23) $\displaystyle\frac{d^{2}\tilde{a}}{dT^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{\Omega_{+}^{\frac{1}{2}}}\frac{d}{dt}\left[\frac{1}{\Omega_{+}^{\frac{1}{2}}}\frac{d}{dt}\left(a\Omega_{-}^{\frac{1}{2}}\right)\right]\
,$ (24)

it follows that the combination

$3{R^{i}}_{i}-{R^{t}}_{t}\equiv\frac{6}{\tilde{a}^{2}}\left(\frac{d\tilde{a}}{dT}\right)^{2}\
,$ (25)

can be written as

$3{R^{i}}_{i}-{R^{t}}_{t}\equiv\frac{6}{\Omega_{-}\Omega_{+}}\left[\frac{1}{a}\frac{d}{dt}\left(a\Omega_{-}^{\frac{1}{2}}\right)\right]^{2}\
,$ (26)

and a little algebra leads to

$\displaystyle\frac{1}{a}\frac{d}{dt}\left(a\Omega_{-}^{\frac{1}{2}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\Omega_{-}^{3/2}}\left[\frac{\dot{a}}{a}\left({\tilde{\lambda}^{2}+2\Phi^{2}}\right)+{\left(\tilde{\lambda}+\Phi\right)}\tfrac{\Phi
V_{\phi}}{\dot{\phi}}\right]\ ,$

which allows us to find an expression for $H\equiv\dot{a}/a$ as

$H=\frac{1}{\left({\tilde{\lambda}^{2}+2\Phi^{2}}\right)}\left[-{\left(\tilde{\lambda}+\Phi\right)}\tfrac{\Phi
V_{\phi}}{\dot{\phi}}\pm\Omega_{-}^{2}\sqrt{\frac{\Omega_{+}}{3\epsilon}\left(|\Omega|^{\frac{1}{2}}-\left(\tilde{\lambda}-2\Phi\right)\right)}\right].$
(27)

The square of this quantity is the generalized FLRW equation for EiBI gravity
coupled to a scalar field with a potential $V(\phi)$. The $\pm$ signs in front
of the square root yield expanding (+) or contracting (-) universes and must
be chosen on physical grounds. In the limit $\epsilon\to 0$ the above
expression yields

$H=\pm\sqrt{\frac{\kappa^{2}(\dot{\phi}^{2}+2V)}{6}}\ ,$ (28)

which recovers the FLRW dynamics of GR coupled to a scalar field.

Note that, from Eq.(12), for a scalar field the energy density $\rho_{\phi}$
and pressure $P_{\phi}$ are related to the field variables by

$\displaystyle\rho_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{\dot{\phi}^{2}}{2}+V$ (29) $\displaystyle P_{\phi}$
$\displaystyle\equiv$ $\displaystyle\frac{\dot{\phi}^{2}}{2}-V\ ,$ (30)

which allows to write the quantities inside the functions $\Omega_{\pm}$ of
Eqs. (14) and (15) in terms of $\rho_{\phi}$ and $P_{\phi}$ as

$\displaystyle\tilde{\lambda}+\Phi$ $\displaystyle=$
$\displaystyle\lambda+\epsilon\kappa^{2}\rho_{\phi}$ (31)
$\displaystyle\tilde{\lambda}-\Phi$ $\displaystyle=$
$\displaystyle\lambda-\epsilon\kappa^{2}P_{\phi}\ .$ (32)

##  3 Massless scalar fields

To explicitly solve the gravitational and scalar field equations (18) and (9)
we need to specify a form of the scalar potential $V(\phi)$. For the sake of
simplicity, and to make contact with similar settings in black hole scenarios
Afonso et al. (2017), let us consider first the case of a free scalar field,
$V=0$, for which the scalar field equation (20) has a first integral

$\dot{\phi}=\frac{\dot{\phi}_{0}}{a^{3}(t)}\ ,$ (33)

with $\dot{\phi}_{0}$ an integration constant. This equation is formally the
same as in the GR case, though the scale factor implicitly contains the
$\epsilon$-corrections via the resolution of Eq.(27), smoothly reducing to
their GR values in the limit $\epsilon\to 0$, as follows from Eq.(28).

Since in this case $P_{\phi}=\rho_{\phi}=\dot{\phi}^{2}/2$, the components of
the deformation matrix (13) assume the relatively simple form

$\Omega_{-}=\left(\lambda^{2}-\frac{\rho_{\phi}^{2}}{\rho_{\epsilon}^{2}}\right)^{\frac{1}{2}}\hskip
2.84544pt;\hskip
2.84544pt\Omega_{+}=\frac{\left(\lambda-s\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{3/2}}{\left(\lambda+s\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{1/2}}\
,$ (34)

where we have defined the critical density
$\rho_{\epsilon}=1/(\kappa^{2}|\epsilon|)$ and $s=\pm 1$ denotes the sign of
$\epsilon$. The corresponding Hubble function takes the form

$H^{2}=\frac{\left(\lambda^{2}-\frac{\rho^{2}_{\phi}}{\rho^{2}_{\epsilon}}\right)^{3/2}\left(\lambda-s\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{2}}{3s\epsilon\left(\lambda^{2}+2\frac{\rho^{2}_{\phi}}{\rho^{2}_{\epsilon}}\right)^{2}}\left(\left(\lambda-s\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)\sqrt{\lambda^{2}-\frac{\rho^{2}_{\phi}}{\rho^{2}_{\epsilon}}}-\lambda+2s\frac{\rho_{\phi}}{\rho_{\epsilon}}\right).$
(35)

As mentioned before, at low densities as compared to the scale
$\rho_{\epsilon}$, this equation recovers the GR expression, which can be
solved leading to the approximate solution $a(t)=H_{0}t^{1/3}$. The Hubble
factor here can be conveniently set to
$H_{0}^{2}=3\kappa^{2}\dot{\phi}^{2}_{0}/2$ to adjust the integration constant
$\dot{\phi}_{0}$ in order to reproduce the standard cosmological evolution of
GR coupled to a massless scalar field. At higher densities, however, the
dynamics strongly departs from that of GR and it is evident that when the
energy density of the scalar field approaches its maximum value,
$\rho_{\phi}=\lambda\rho_{\epsilon}$, the Hubble function vanishes. This
corresponds to a minimum value of the expansion factor for both branches of
solutions $s=\pm 1$, but there are two different mechanisms by which the
initial Big Bang singularity is avoided, which we discuss next. For
simplicity, from now on we focus on asymptotically flat configurations,
$\lambda=1$.

Figure 1: Representation of the Hubble function
$\sqrt{|\epsilon|}H(\rho_{\phi})$ vs. the energy density for EiBI gravity
(solid) and GR (dashed) with massless scalar field as a function of the ratio
$\rho_{\phi}/\rho_{\epsilon}$. Note how the trajectories of the EiBI theory
are bounded, while that of GR (dashed black) is open. Bouncing solutions
(blue) reach the maximum density forming a $\pi/2$ angle with the horizontal
axis, while the loitering branch (red) reaches it tangentially.

For the $s=-1$ branch, at the critical density $\rho_{\phi}=\rho_{\epsilon}$,
the Hubble factor in Eq.(35) and its derivative behave as

$\displaystyle H(\rho_{\phi})$ $\displaystyle\approx$
$\displaystyle\pm\frac{2^{7/4}}{3|\epsilon|^{1/2}}\left(1-\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{3/4}$
(36) $\displaystyle\frac{dH}{d\rho_{\phi}}$ $\displaystyle\approx$
$\displaystyle\mp\frac{1}{2^{\frac{1}{4}}|\epsilon|^{\frac{1}{2}}\rho_{\epsilon}\left(1-\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{1/4}}\
.$ (37)

This indicates that at the maximum achievable density,
$\rho_{\phi}=\rho_{\epsilon}$, the Hubble factor vanishes while its derivative
goes to infinity. This implies that the universe contracts to a minimum
(maximum) value of the expansion factor (energy density), before re-expanding.
This the typical behaviour expected in standard bouncing solutions with a
transition from a contracting phase to an expanding one, where the universe
could even undergo a sequence of cyclic cosmologies with both phases. In Fig.
1 we plot the form of $H(\rho)$ (blue curves) according to Eq.(35) to show
that while at late-times (i.e. low densities, $\rho_{\phi}\ll\rho_{\epsilon}$)
these solutions converge to those of GR (in agreement to Eq.(28)), at high
densities ($\rho_{\phi}/\rho_{\epsilon}\to 1$) they depart form the standard
Big Bang singularity of GR.

For the $s=+1$ branch we find that the evolution of the Hubble factor as we
approach the maximum density $\rho_{\phi}=\rho_{\epsilon}$ is given instead by

$\displaystyle H(\rho_{\phi})$ $\displaystyle\approx$
$\displaystyle\pm\frac{2^{3/4}\sqrt{\frac{1}{\epsilon}}\left(1-\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{7/4}}{3\sqrt{3}}$
(38) $\displaystyle\frac{dH}{d\rho_{\phi}}$ $\displaystyle\approx$
$\displaystyle\mp\frac{7\sqrt{\frac{1}{\epsilon}}\left(1-\frac{\rho_{\phi}}{\rho_{\epsilon}}\right)^{3/4}}{6\sqrt[4]{2}\sqrt{3}\rho_{\epsilon}}\
,$ (39)

which means that the Hubble factor vanishes there as well, but instead of
being divergent its derivative takes a finite value (which is actually zero).
This implies that the expansion factor reaches a fixed value $a(t)=a_{m}$,
corresponding to an asymptotically Minkowski past, and starts expanding as we
move forward in time after some reference time $t=t_{p}$. In Fig. 1 we depict
this behaviour (red curves) of these so-called loitering solutions, where the
qualitative differences with the bouncing solutions are manifest. Again, at
late times (low densities) the standard GR evolution is recovered. Therefore
we see that a free massless scalar field coupled to EiBI gravity is able to
yield nonsingular evolutions in both $s=\pm 1$ branches according to these two
mechanisms, something not possible within GR.

##  4 Massive scalar field

Let us now consider a massive scalar field with a potential of the form

$V(\phi)=\frac{1}{2}\mu^{2}\phi^{2}\ ,$ (40)

where $\mu$ is a constant. This choice for the potential, besides being the
simplest form, corresponds to the natural mass associated to a canonical
scalar field. In this case, the Hubble function (27) develops an explicit
dependence on both $\rho_{\phi}$ and $P_{\phi}$ [to lighten the notation from
now on we shall drop the label $\phi$ in these quantities] via Eqs.(31) and
(32). Since the resulting expression does not provide any useful insight, we
will not write it explicitly here. Instead, for the sake of comparison with
the massless case of the previous section, we find it more illustrative to
provide a parametric representation of $H(t)$ versus $\rho(t)$. These
functions are obtained by numerically integrating the second-order equation
for $\ddot{a}$ that follows from (22) together with its corresponding right-
hand side. The result of the integration is used to construct the quantity
$\dot{a}/a$, which is then compared with the formula (27) to check the
consistency of the numerical integration. Again we split our discussion of the
corresponding results into the $s=\pm 1$ cases.

Figure 2: Parametric plot of $\sqrt{|\epsilon|}H(t)$ as a function of
$\rho(t)/\rho_{\epsilon}$ for small masses when $s=-1$ (bouncing solutions).
The (reduced) mass parameter is taken such that
$\sigma\equiv\mu^{2}/\rho_{\epsilon}$.  Figure 3: Parametric plot of
$\sqrt{|\epsilon|}H(t)$ as a function of $\rho(t)$ for larger values of the
reduced mass $\sigma\equiv\mu^{2}/\rho_{\epsilon}$ when $s=-1$. As compared to
Fig. 2, in this case the development of fish-like structures is clearly
visible as $\sigma$ grows large enough. Figure 4: Expansion factor $a(t)$
representing bouncing solutions for various values of the reduced mass
parameter $\sigma\equiv\mu^{2}/\rho_{\epsilon}$ when $s=-1$. Figure 5: Time
derivative of the expansion factor representing bouncing solutions for various
values of the reduced mass parameter $\sigma\equiv\mu^{2}/\rho_{\epsilon}$.
The colors of the curves are the same as in Fig. 4. Note that the blue and
magenta curves will bounce at some point and begin a growing oscillatory
trajectory similar to that of the red dashed curve, being this a generic
behavior of all these solutions.

In Figs. 2 and 3, we show the function $H(\rho(t))$ from Eq.(27) with the help
of (14) and (15) particularized to the potential above, for solutions
corresponding to several values of the reduced mass
$\sigma\equiv\mu^{2}/\rho_{\epsilon}$ in the case $s=-1$. The corresponding
expansion factors $a(t)$ [obtained from integrating the previous function $H$]
appear in Fig. 4, while in Fig. 5 we show the function $\dot{a}(t)$ of those
same examples and in Fig. 6 their Hubble functions. As one can see from all
these plots, the case $s=-1$ still represents bouncing solutions for all the
mass range explored, which goes from $\sigma\approx 10^{-8}$ up to
$\sigma\approx 0.4$, with little variations with respect to the massless
scenario for masses as high as $\sigma\approx 10^{-3}$ (see Fig. 2). For
higher masses, the egg-shaped Hubble function develops a fish-like structure,
with asymmetric fins. This asymmetry is also manifest in the contracting
branch of the expansion factor (see Fig. 4), which becomes increasingly
asymmetric as the mass parameter grows from zero. A similar behaviour can be
observed in the density profiles of the scalar field depicted in Fig. 7.

Figure 6: Hubble function $H\equiv\dot{a}/a$ for various values of the reduced
mass parameter $\sigma\equiv\mu^{2}/\rho_{\epsilon}$ of the $s=-1$ case.
Figure 7: Energy density (solid curves) as a function of time superimposed
with its corresponding expansion factor $a(t)$ (dashed curves) for some
bouncing solutions of the $s=-1$ case. Figure 8: Parametric plot of $H(t)$ as
a function of $\rho(t)$ for small values of the reduced mass $\sigma$ of the
$s=+1$ case. Figure 9: Expansion factor representing new bouncing solutions
for various values of the reduced mass parameter $\sigma$ when $s=+1$.  Figure
10: Time derivative of the expansion factor representing bouncing solutions
for the reduced mass parameter $\sigma$ when $s=+1$. Figure 11: Hubble
function $\dot{a}/a$ for various values of the reduced mass parameter $\sigma$
when $s=+1$. Figure 12: Energy density (solid curves) as a function of time
superimposed with its corresponding expansion factor $a(t)$ (dashed curves).
Figure 13: Zoom in of the energy density (solid curves) as a function of time
superimposed with its corresponding expansion factor $a(t)$ (dashed curves).

Focusing now on the case with $s=+1$, which represents loitering solutions in
the massless case, Fig. 8 shows that the function $H(\rho(t))$ is now
qualitatively different (compare to Fig. 1). Indeed, when $\sigma=0$, $H(t)$
never crosses the horizontal axis as one approaches the maximum allowed energy
density. However, for $\sigma\neq 0$ one either has an expanding branch which
starts from an almost constant $a(t)=a_{m}$ initial phase or a contracting
phase that ends up in an almost constant final $a(t)=a_{m}$. Now we observe
that all solutions develop a bounce at some high density, effectively crossing
the axis and establishing a continuous transition from a contracting phase to
an expanding one. The loitering phase, with an almost constant $a(t)$, may
last for a long period after an initial contraction but it will always end up
in an expanding branch. This is illustrated in Fig. 9, where the blue dashed
curve represents a phase with $a(t)$ almost constant for a long time (very low
mass). The orange curve has essentially the same future behaviour as the blue
curve but its loitering phase (almost constant expansion factor) does not last
so much, exhibiting a previously contracting phase that arises abruptly.
Similar behaviours arise for larger masses, and one verifies that the
instabilities in the expansion factor always lead to a bounce, never finding
fully collapsed solutions. This is a remarkable property of the EiBI model,
since it always avoids the development of singular solutions.

For completeness, the behaviour of $\dot{a}(t)$ and $H(t)$ for this case
$s=+1$ are shown in Figs. 10 and 11. The corresponding energy density profiles
appear in Figs. 12 and 13. It is amusing to see how for low mass
configurations (orange solid curve, for instance) the energy density at early
times is very low (large universe with very diluted energy) until it rapidly
increases to reach a maximum, where it stays for some time at an almost
constant value (loitering phase) with a slight decay (decompression) as we
move forward in time. Then the density suddenly drops again giving rise to an
expanding phase. For larger values of the scalar field mass, this process can
be significantly deformed but the qualitative features remain.

##  5 Observational Constraints

In this section we shall put to observational test the EiBI cosmologies
considered in the previous sections by using the latest observational data to
constrain the EiBI parameter $\epsilon$ preventing the singularity in the
early Universe.

###  5.1 BBN

In the first sections we began with the original action with the scalar field
inside, and shifted the solution to other frame. This frame, as we showed,
removes the singularity of the solution and for $\epsilon\rightarrow 0$ it
recovers GR. Here, we put the density of the $\Lambda$CDM that produces in the
"Einstein frame" a solution that still removes the singularity and recovers GR
for $\epsilon\rightarrow 0$. Doing the same procedure for a theory with matter
fields $\rho_{m}$ in the original frame gives the density:

$\tilde{\rho}=\frac{8(3p\epsilon-1)^{2}(3\rho\epsilon+1)\left(2\sqrt{1-3p\epsilon}\sqrt{3\rho\epsilon+1}-3\epsilon\left(p\left(2\sqrt{1-3p\epsilon}\sqrt{3\rho\epsilon+1}-3\right)-\rho\right)-2\right)}{3\epsilon\left((3\rho\epsilon+1)\left(3\epsilon\left(ap^{\prime}+4p\right)-4\right)+3a\epsilon(3p\epsilon-1)\rho^{\prime}\right)^{2}}$
(41)

in the physical frame. For $\epsilon\rightarrow 0$ the relation gives
$\tilde{\rho}=\rho$, as expected.

In order to constrain the EiBI parameter $\epsilon$ we use different
measurements from our Universe. The strongest one is from the Big Bang
Nucleosynthesis (BBN) where $z\sim 10^{9}$. In principle the condition for the
BBN constraint is Solà et al. (2017); Barrow et al. (2021):

$\left|\frac{\rho_{EiBI}-\rho_{\Lambda CDM}}{\rho_{\Lambda CDM}}\right|<10\%$
(42)

In this era the matter dominance is radiation, so we take the density and the
pressure to be:

$\rho=\frac{\Omega_{r}}{a^{4}},\quad p=\frac{\rho}{3}.$ (43)

The density for the EiBI theory gives:

$\tilde{\Omega}=\frac{\left(a^{4}-\Omega_{r}\epsilon\right)^{2}\left(a^{4}+3\Omega_{r}\epsilon\right)\left(a^{4}\left(\sqrt{1-\frac{\Omega_{r}\epsilon}{a^{4}}}\sqrt{\frac{3\Omega_{r}\epsilon}{a^{4}}+1}-1\right)+\Omega_{r}\epsilon\left(3-\sqrt{1-\frac{\Omega_{r}\epsilon}{a^{4}}}\sqrt{\frac{3\Omega_{r}\epsilon}{a^{4}}+1}\right)\right)}{3\epsilon\left(a^{8}+3\Omega_{r}^{2}\epsilon^{2}\right)^{2}}.$
(44)

From the condition Eq.(42) we get the bound444Recall that restoring dimensions
we have $\epsilon\to\kappa^{2}\epsilon$, which represents the inverse of a
matter density. Accordingly, the bound on $\epsilon$, which represents a
squared length, can be written also as $\left|\epsilon\right|\lesssim
5.53\cdot 10^{-8}\text{ m}^{2}$. :

$\left|\kappa^{2}\epsilon\right|\lesssim 1.03\cdot 10^{-33}.$ (45)

where we took $\Omega_{r}\sim 10^{-4}$. It is important to note that this
bound improves any previous existing bounds on $\epsilon$ by several orders of
magnitude. In fact, the most stringent constraints so far came from elementary
particles scattering experiments Jiménez et al. (2021); Latorre et al. (2018);
Delhom et al. (2020) and implied $\kappa^{2}\epsilon\lesssim 1.86\cdot
10^{-28}$ (or equivalently, $\epsilon\lesssim 10^{-2}\text{ m}^{2}$).

We now proceed to describe the observational data sets along with the relevant
statistics in constraining the model, while we use the first constraint from
the BBN above. The matter fields considered are the dark energy
$\Omega_{\Lambda}$, dark matter $\Omega_{m}$, and radiation $\Omega_{r}$
components. The corresponding energy density reads

$\frac{\rho}{\rho_{c}}=\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}}+\Omega_{\Lambda}\
.$ (46)

Inserting this equation in the expression of the Hubble factor (27), and
assuming a small EiBI parameter, $|R_{\mu\nu}|\ll\epsilon^{-1}$, yields the
extended Friedman equation:

$\frac{H^{2}}{H_{0}^{2}}=\frac{\Omega_{m}}{a^{3}}+\frac{\Omega_{r}}{a^{4}}+\Omega_{\Lambda}+\epsilon\left(\frac{9\Omega_{m}^{2}}{8a^{6}}+\frac{\Omega_{m}\Omega_{r}}{a^{7}}-\frac{2\Omega_{\Lambda}\Omega_{r}}{a^{4}}\right)+\mathcal{O}(\epsilon^{2}),$
(47)

where the EiBI corrections to GR solutions are apparent.

###  5.2 Direct measurements of the Hubble expansion

Cosmic Chronometers (CC): This data set exploits the evolution of differential
ages of passive galaxies at different redshifts to directly constrain the
Hubble parameter Jimenez and Loeb (2002). We use uncorrelated 30 CC
measurements of $H(z)$ discussed in Moresco et al. (2012a, b); Moresco (2015);
Moresco et al. (2016). The corresponding $\chi^{2}_{H}$ function reads:

$\chi^{2}_{H}=\sum_{i=1}^{30}\left(\frac{H_{i}-H_{pred}(z_{i})}{\Delta
H_{i}}\right)^{2},$ (48)

where $H_{i}$ is the observed Hubble rates at redshift $z_{i}$ ($i=1,...,N$)
and $H_{pred}$ is the predicted one from the model.

###  5.3 Standard Candles

As Standard Candles (SC) we use measurements of the Pantheon type Ia supernova
Scolnic et al. (2018), that were collected in Anagnostopoulos et al. (2020),
as well as quasars Roberts et al. (2017) and gamma ray bursts Demianski et al.
(2017). The model parameters are fitted by comparing the observed
$\mu_{i}^{obs}$ value to the theoretical $\mu_{i}^{th}$ value of the distance
moduli, which are the logarithms:

$\mu=m-M=5\log_{10}(D_{L})+\mu_{0}\ ,$ (49)

where $m$ and $M$ are the apparent and absolute magnitudes and
$\mu_{0}=5\log\left(H_{0}^{-1}/Mpc\right)+25$ is the nuisance parameter that
has been marginalized. The luminosity distance is defined by

$\displaystyle D_{L}(z)$ $\displaystyle=$
$\displaystyle\frac{c}{H_{0}}(1+z)\int_{0}^{z}\frac{dz^{\prime}}{E(z^{\prime})}\
,$ (50)

(here $\Omega_{k}=0$, i.e., a flat space-time and $E(z)$ is the dimensionless
Hubble parameter). For the SnIa data the covariance matrix is not diagonal and
the distance modulus is given as $\mu_{i}=\mu_{B,i}-\mathcal{M}$, where
$\mu_{B,i}$ is the apparent magnitude at maximum in the rest frame for
redshift $z_{i}$ and $\mathcal{M}$ is treated as a universal free parameter,
quantifying various observational uncertainties Scolnic et al. (2018).
Following standard lines, the chi-square function of the standard candles is
given by

$\chi^{2}_{\text{SC}}\left(\phi^{\nu}_{\text{s}}\right)={\mu}_{\text{s}}\,{C}_{\text{s},\text{cov}}^{-1}\,{\mu}_{\text{s}}^{T}\,,$
(51)

where
${\mu}_{\text{s}}=\\{\mu_{1}-\mu_{\text{th}}(z_{1},\phi^{\nu})\,,\,...\,,\,\mu_{N}-\mu_{\text{th}}(z_{N},\phi^{\nu})\\}$
and the subscript ‘s’ denotes SnIa and QSOs.

###  5.4 Baryon acoustic oscillations

We use uncorrelated data points from different Baryon Acoustic Oscillations
(BAO). BAO are a direct consequence of the strong coupling between photons and
baryons in the pre-recombination epoch. After the decoupling of photons, the
over densities in the baryon fluid evolved and attracted more matter, leaving
an imprint in the two-point correlation function of matter fluctuations with a
characteristic scale of around $r_{d}\approx 147$Mpc that can be used as a
standard ruler and to constrain cosmological models. Studies of the BAO
feature in the transverse direction provide a measurement of
$D_{H}(z)/r_{d}=c/H(z)r_{d}$, with the comoving angular diameter distance Hogg
et al. (2020); Martinelli et al. (2020):

$D_{M}=\int_{0}^{z}\frac{c\,dz^{\prime}}{H(z^{\prime})}\ .$ (52)

The angular diameter distance $D_{A}=D_{M}/(1+z)$ and the quantity
$D_{V}(z)/r_{d}$ with

$D_{V}(z)\equiv[zD_{H}(z)D_{M}^{2}(z)]^{1/3}$ (53)

are a combination of the BAO peak coordinates above. The surveys provide the
values of the measurements at some effective redshift. We employ the BAO data
points collected in Benisty and Staicova (2020) from Percival et al. (2010);
Beutler et al. (2011); Busca et al. (2013); Anderson et al. (2013); Seo et al.
(2012); Ross et al. (2015); Tojeiro et al. (2014); Bautista et al. (2018); de
Carvalho et al. (2018); Ata et al. (2018); Abbott et al. (2019); Molavi and
Khodam-Mohammadi (2019), in the redshift range $0.106<z<2.34$. The BAO scale
is set by the sound horizon at the drag epoch $z_{d}\approx 1060$ when photons
and baryons decouple. In our analysis we used $r_{d}$ as independent
parameter.

The BAO data represent the absolute distance measurements in the Universe.
From the measurements of correlation function or power spectrum of large scale
structure, we can use the BAO signal to estimate the distance scales at
different redshifts. In practice, the BAO data are analyzed based on a
fiducial cosmology and the sound horizon at drag epoch. We calculate the
redshift of the drag epoch as:

$z_{d}=\frac{1291\omega_{m}^{0.251}}{1+0.659\omega_{m}^{0.828}}[1+b_{1}\omega_{b}^{b_{2}}],$
(54)

where

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle
0.313\omega_{m}^{-0.419}[1+0.607\omega_{m}^{0.674}],$ (55) $\displaystyle
b_{2}$ $\displaystyle=$ $\displaystyle 0.238\omega_{m}^{0.223}.$ (56)

The uncorrelation of this data set yields the corresponding $\chi^{2}$:

$\chi^{2}_{BAO}=\sum_{i=1}^{17}\left(\frac{D_{i}-D_{pred}(z_{i})}{\Delta
D_{i}}\right)^{2},$ (57)

where $D_{i}$ is the observed distant module rates at redshift $z_{i}$
($i=1,...,N$) and $D_{pred}$ is the predicted one from the model.

###  5.5 CMB distant priors

We use the CMB distant priors that based on the latest CMB measurements
(Aghanim et al., 2020). Its contribution in the likelihood analysis is
expressed in terms of the compressed form with CMB shift parameters:

$\displaystyle R$ $\displaystyle\equiv$
$\displaystyle\sqrt{\Omega_{m}H_{0}^{2}}r(z_{*})/c,$ (58) $\displaystyle
l_{a}$ $\displaystyle\equiv$ $\displaystyle\pi r(z_{*})/r_{s}(z_{*}),$ (59)

where $r_{s}(z)$ is the comoving sound horizon at redshift $z$, and $z_{*}$ is
the redshift to the photon-decoupling surface. These two CMB shift parameters
together with $\omega_{b}=\Omega_{b}h^{2}$ and spectral index of the
primordial power spectrum $n_{s}$ can give an efficient summary of CMB data
for the dark energy constraints.

The comoving sound horizon is given by

$r_{s}(z)=\frac{c}{H_{0}}\int_{0}^{a}\frac{da^{\prime}}{\sqrt{3(1+\bar{R_{b}}a^{\prime})a^{\prime
4}E^{2}(z^{\prime})}}.$ (60)

The radiation term in the expression of $E(z)$ for the CMB data analysis
should not be ignored. It can be determined by the matter-radiation equality
relation $\Omega_{r}=\Omega_{m}/(1+z_{\text{eq}})$, and
$z_{\text{eq}}=2.5\times 10^{4}\omega_{m}(T_{\text{CMB}}/2.7\text{K})^{-4}$,
where $\omega_{m}=\Omega_{m}h^{2}$. The sound speed is
$c_{s}=1/\sqrt{3(1+\bar{R_{b}}a)}$, with $\bar{R_{b}}a=3\rho_{b}/(4\rho_{r})$,
and $\bar{R_{b}}=31500w_{b}(T_{\text{CMB}}/2.7\text{K})^{-4}$. We assume the
CMB temperature $T_{\text{CMB}}=2.7255\text{K}$. The redshift $z_{*}$ can be
calculated by the fitting formula:

$z_{*}=1048[1+0.00124\omega_{b}^{-0.738}][1+g_{1}\omega_{m}^{g_{2}}],$ (61)

where

$g_{1}=\frac{0.0783\omega_{b}^{-0.238}}{1+39.5\omega_{b}^{0.763}},\qquad
g_{2}=\frac{0.560}{1+21.1\omega_{b}^{1.81}}.$ (62)

###  5.6 Direct Detection of the Hubble parameter

We include the latest measurement of the Hubble parameter:

$H_{0}=(73.2\pm 1.3)\text{km/s/Mpc}$ (63)

reported by Riess et al. (2021). The measurement presents an expanded sample
of 75 Milky Way Cepheids with Hubble Space Telescope (HST) photometry and Gaia
EDR3 parallaxes which uses the extragalactic distance ladder in order to
recalibrate and refine the determination of the Hubble constant. The
combination is related via the relation

$\chi^{2}_{Hub}=\left(\frac{H_{0}-73.2}{1.3}\right)^{2}.$ (64)

The $\chi^{2}_{Hub}$ estimates the deviation from the latest measurement of
the Hubble constant.

###  5.7 Joint analysis and model selection

In order to perform a joint statistical analysis of $4$ cosmological probes we
need to use the total likelihood function, consequently the
$\chi^{2}_{\text{tot}}$ expression is given by

$\chi_{\text{tot}}^{2}=\chi_{BBN}^{2}+\chi_{CMB}^{2}+\chi_{H}^{2}+\chi_{SC}^{2}+\chi_{BAO}^{2}+\chi^{2}_{Hub}.$
(65)

In order to perform a joint statistical analysis of these four cosmological
probes we need to use the total likelihood function. Regarding the problem of
likelihood maximization, we use an affine-invariant Markov Chain Monte Carlo
sampler Foreman-Mackey et al. (2013), as it is implemented within the open-
source packaged $Polychord$ Handley et al. (2015) with the $GetDist$ package
Lewis (2019) to present the results. The prior we choose is with a uniform
distribution, where $\Omega_{m}\in[0.;1.]$,
$\Omega_{\Lambda}\in[0.;1-\Omega_{m}]$, $H_{0}\in[50;100]$Km/sec/Mpc. For the
EiBI parameter $\epsilon$ we set the range
$|\kappa^{2}\epsilon|\in[0.;10^{-33}]$.

Model | $H_{0}(Km/sec/Mpc)$ | $\Omega_{m}$ | $\Omega_{\Lambda}$ | $\epsilon(10^{-8})$ | $\Delta B_{ij}$
---|---|---|---|---|---
EiBI | $64.99\pm 0.556$ | $0.291\pm 0.0074$ | $0.689\pm 0.0084$ | $\left(5.76\pm 3.53\right)$ | $0.51$
$\Lambda$CDM | $64.94\pm 0.525$ | $0.291\pm 0.008$ | $0.6891\pm 0.00865$ | - | -

Figure 14: The posterior distribution for the simplest case of the EiBI and
for $\Lambda$CDM model (blue curve). The ratios of the matter density
$\Omega_{m}$, dark energy $\Omega_{\Lambda}$. The final results for
cosmological parameters for the EiBI and the $\Lambda$CDM models are
summarized in the table. In order to compare the models we calculate the Bayes
factor.

The posterior distribution of $H_{0}$ vs. the EiBI parameter $\epsilon$ is
presented in Fig. 14 and 15. The Hubble parameter is $H_{0}=64.99\pm
0.556$km/sec/Mpc which is in between the Planck estimation of the Hubble
function Aghanim et al. (2020) and the latest SH0ES measurement Riess et al.
(2021). The dark matter component is $\Omega_{m}=0.291\pm 0.0074$ and the dark
energy component is being $\Omega_{\Lambda}=0.689\pm 0.0084$. The fit for the
EiBI parameter gives $\kappa^{2}\epsilon=\left(4.764\pm 3.074\right)\cdot
10^{-33}$. We point out that the value of $\epsilon=0$ (which corresponds to
GR) is very close to the mean value. Therefore, for
$|\epsilon|\lesssim\left(5.76\pm 3.53\right)\cdot 10^{-8}$m2 the theory is
able to fit to the data while being able to successfully remove the Big Bang
singularity. We point out again that this improves by several orders of
magnitude the previous strongest constraints for EiBI gravity’s parameter as
obtained from particle physics experiments Jiménez et al. (2021); Latorre et
al. (2018); Delhom et al. (2020).

We finally point out that from the Bayes factor difference between the models
that reads $\Delta B_{ij}=0.51$ we get an indistinguishable difference for the
$\Lambda$CDM model. Therefore, EiBI gravity with the constraint above for its
parameter when coupled to a scalar yields very similar properties for the
$\Lambda$CDM model, but prevents the initial singularity.

Figure 15: The posterior distribution for the Hubble parameter $H_{0}$ vs. the
EiBI parameter $\epsilon$.

##  6 Conclusion and discussion

In this paper we have analyzed homogeneous and isotropic cosmological
solutions in the context of Eddington-inspired Born-Infeld gravity coupled to
a single scalar field. In the massless case we discussed the existence and
properties of such solutions by direct resolution of the modified
gravitational field equations, which depend on the sign of the EiBI parameter.
For the bouncing solutions ($s=-1$) one finds that at the maximum achievable
density of the scalar field the expansion factor gets to a minimum non-
vanishing value (vanishing Hubble factor) undergoing a transition from a
contracting phase to an expanding one (or vice versa). In the loitering
solutions ($s=+1$) the cosmological evolution starts from an asymptotically
Minkowski past with a minimum constant value of the expansion factor, such
that at a given time a canonical cosmological expansion phase is triggered.
These results are consistent with the expectancy raised in the seminal paper
Banados and Ferreira (2010) and the numerical evidence found in more recent
works Beltran Jimenez et al. (2017) when considering perfect fluids as the
matter source.

Though the massless case was known to have only two types of solutions,
namely, bouncing and loitering, some kind of instability was expected for the
loitering case beyond the massless scenario. Until now it was not known if
they would lead to singular solutions of if they would remain nonsingular. Our
analysis puts forward that the model is absolutely robust against
singularities and that the loitering case gives rise to a new type of
solutions in which an unstable Minkowskian phase (almost constant expansion
factor) is possible in between a contraction phase and an expansion phase. The
duration of this Minkowskian period (see Fig. 9) depends on model parameters
and the initial conditions. No singular solutions are found in our analysis,
neither in the massless nor in the massive case. All the nonsingular solutions
recover the late-time GR cosmological evolution.

Getting into the details, for the $s=-1$ (bouncing) branch we found that for
low masses the solutions closely resemble the behaviour of those of the
massless case. However, for masses above a certain threshold the parametric
representation of the Hubble function versus the energy density develops a
kind of fish-like asymmetric fins at certain time in the cosmological
evolution, a feature not present in the massless case (compare Fig. 1 with
Figs. 2 and 3). The $s=+1$ (loitering) branch also exhibits unusual features,
showing that during the loitering (almost Minkowskian) phase the energy
density can grow significantly while the Hubble function remains close to zero
and positive (see Fig. 8). This represents a kind of intermediate state in a
continuous transition from a contracting phase to a final expansion one. This
behaviour is triggered by the energy density stored in the massive scalar
field, whose drop marks the end of the loitering phase and the beginning of
the expansion phase (see Figs. 12 and 13).

To conclude, the results obtained in this paper further support the
effectiveness of EiBI gravity to remove singularities in cosmological and
astrophysical scenarios when coupled to different matter sources. In the
present scenario with scalar fields this can be done while being consistent
with the latest cosmological observations provided that the EiBI parameter is
kept roughly within $\left|\epsilon\right|\lesssim 5\cdot 10^{-8}\text{
m}^{2}$. One important aspect of the obtained solutions not tackled here is
whether they are stable under tensorial perturbations or not. We point out
that in the present case the equation of state is of the form
$\omega=\omega(t)$ since the energy density and the pressure of the scalar
field (in the massive case) are not trivially related. In turn this should
have an impact on the behaviour of the quantities prone to the development of
instabilities in the tensorial perturbations equation, such as the sound speed
Beltran Jimenez et al. (2017); Escamilla-Rivera et al. (2012). Therefore, one
would need to carry out such an analysis in order to further support the
feasibility of this theory and their associated solutions to replace the Big
Bang singularity by observationally viable nonsingular cosmologies.

###### Acknowledgements.

D.B. gratefully acknowledges the support from the Blavatnik and the Rothschild
fellowship. G.J.O. is no longer funded by the Ramon y Cajal contract
RYC-2013-13019 (Spain) because he became permanent. D.R.-G. is funded by the
_Atracción de Talento Investigador_ programme of the Comunidad de Madrid
(Spain) No. 2018-T1/TIC-10431, and acknowledges further support from the
Ministerio de Ciencia, Innovación y Universidades (Spain) project No.
PID2019-108485GB-I00/AEI/10.13039/501100011033, and the FCT projects No.
PTDC/FIS-PAR/31938/2017 and PTDC/FIS-OUT/29048/2017. This work is supported by
the Spanish Grant FIS2017-84440-C2- 1-P funded by MCIN/AEI/
10.13039/501100011033 “ERDF A way of making Europe”, Grant
PID2020-116567GB-C21 funded by MCIN/AEI/10.13039/501100011033, the project
PROMETEO/2020/079 (Generalitat Valenciana), the project i-COOPB20462 (CSIC)
and the Edital 006/2018 PRONEX (FAPESQ-PB/CNPQ, Brazil, Grant 0015/2019). This
article is based on work from COST Action CA18108, supported by COST (European
Cooperation in Science and Technology

## References

  * Amendola et al. (2013) Amendola, L.; others.  Cosmology and fundamental physics with the Euclid satellite.  Living Rev. Rel. 2013, 16, 6, [arXiv:astro-ph.CO/1206.1225].  doi:black10.12942/lrr-2013-6.
  * Aghanim et al. (2020) Aghanim, N.; others.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys. 2020, 641, A6, [arXiv:astro-ph.CO/1807.06209].  doi:black10.1051/0004-6361/201833910.
  * Starobinsky (1979) Starobinsky, A.A.  Spectrum of relict gravitational radiation and the early state of the universe.  JETP Lett. 1979, 30, 682–685.  [,767(1979)].
  * Starobinsky (1980) Starobinsky, A.A.  A New Type of Isotropic Cosmological Models Without Singularity.  Phys. Lett. 1980, 91B, 99–102.  771 (1980), doi:black10.1016/0370-2693(80)90670-X.
  * Guth (1981) Guth, A.H.  The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems.  Phys. Rev. 1981, D23, 347–356.  [Adv. Ser. Astrophys. Cosmol.3,139(1987)], doi:black10.1103/PhysRevD.23.347.
  * Albrecht and Steinhardt (1982) Albrecht, A.; Steinhardt, P.J.  Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking.  Phys. Rev. Lett. 1982, 48, 1220–1223.  [Adv. Ser. Astrophys. Cosmol.3,158(1987)], doi:black10.1103/PhysRevLett.48.1220.
  * Mukhanov and Chibisov (1981) Mukhanov, V.F.; Chibisov, G.V.  Quantum Fluctuations and a Nonsingular Universe.  JETP Lett. 1981, 33, 532–535.  [Pisma Zh. Eksp. Teor. Fiz.33,549(1981)].
  * Guth and Pi (1982) Guth, A.H.; Pi, S.Y.  Fluctuations in the New Inflationary Universe.  Phys. Rev. Lett. 1982, 49, 1110–1113.  doi:black10.1103/PhysRevLett.49.1110.
  * Linde (1982) Linde, A.D.  A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems.  Phys. Lett. 1982, 108B, 389–393.  [Adv. Ser. Astrophys. Cosmol.3,149(1987)], doi:black10.1016/0370-2693(82)91219-9.
  * Barrow and Cotsakis (1988) Barrow, J.D.; Cotsakis, S.  Inflation and the Conformal Structure of Higher Order Gravity Theories.  Phys. Lett. 1988, B214, 515–518.  doi:black10.1016/0370-2693(88)90110-4.
  * Barrow (1988) Barrow, J.D.  The Premature Recollapse Problem in Closed Inflationary Universes.  Nucl. Phys. B 1988, 296, 697–709.  doi:black10.1016/0550-3213(88)90040-5.
  * Elizalde et al. (2008) Elizalde, E.; Nojiri, S.; Odintsov, S.D.; Saez-Gomez, D.; Faraoni, V.  Reconstructing the universe history, from inflation to acceleration, with phantom and canonical scalar fields.  Phys. Rev. D 2008, 77, 106005, [arXiv:hep-th/0803.1311].  doi:black10.1103/PhysRevD.77.106005.
  * Ratra and Peebles (1988) Ratra, B.; Peebles, P.J.E.  Cosmological Consequences of a Rolling Homogeneous Scalar Field.  Phys. Rev. 1988, D37, 3406.  doi:black10.1103/PhysRevD.37.3406.
  * Caldwell et al. (1998) Caldwell, R.R.; Dave, R.; Steinhardt, P.J.  Cosmological imprint of an energy component with general equation of state.  Phys. Rev. Lett. 1998, 80, 1582–1585, [arXiv:astro-ph/astro-ph/9708069].  doi:black10.1103/PhysRevLett.80.1582.
  * Kehayias and Scherrer (2019) Kehayias, J.; Scherrer, R.J.  New generic evolution for $k$ -essence dark energy with $w\approx-1$.  Phys. Rev. 2019, D100, 023525, [arXiv:gr-qc/1905.05628].  doi:black10.1103/PhysRevD.100.023525.
  * Oikonomou and Chatzarakis (2020) Oikonomou, V.K.; Chatzarakis, N.  The Phase Space of $k$-Essence $f(R)$ Gravity Theory.  Nucl. Phys. 2020, B956, 115023, [arXiv:gr-qc/1905.01904].  doi:black10.1016/j.nuclphysb.2020.115023.
  * Chakraborty et al. (2019) Chakraborty, A.; Ghosh, A.; Banerjee, N.  Dynamical systems analysis of a k -essence model.  Phys. Rev. 2019, D99, 103513, [arXiv:gr-qc/1904.10149].  doi:black10.1103/PhysRevD.99.103513.
  * Babichev et al. (2018) Babichev, E.; Ramazanov, S.; Vikman, A.  Recovering $P(X)$ from a canonical complex field 2018.  [arXiv:gr-qc/1807.10281].  [JCAP1811,023(2018)], doi:black10.1088/1475-7516/2018/11/023.
  * Zlatev et al. (1999) Zlatev, I.; Wang, L.M.; Steinhardt, P.J.  Quintessence, cosmic coincidence, and the cosmological constant.  Phys. Rev. Lett. 1999, 82, 896–899, [arXiv:astro-ph/astro-ph/9807002].  doi:black10.1103/PhysRevLett.82.896.
  * Caldwell (2002) Caldwell, R.R.  A Phantom menace?  Phys. Lett. 2002, B545, 23–29, [arXiv:astro-ph/astro-ph/9908168].  doi:black10.1016/S0370-2693(02)02589-3.
  * Chiba et al. (2000) Chiba, T.; Okabe, T.; Yamaguchi, M.  Kinetically driven quintessence.  Phys. Rev. 2000, D62, 023511, [arXiv:astro-ph/astro-ph/9912463].  doi:black10.1103/PhysRevD.62.023511.
  * Bento et al. (2002) Bento, M.C.; Bertolami, O.; Sen, A.A.  Generalized Chaplygin gas, accelerated expansion and dark energy matter unification.  Phys. Rev. 2002, D66, 043507, [arXiv:gr-qc/gr-qc/0202064].  doi:black10.1103/PhysRevD.66.043507.
  * Tsujikawa (2013) Tsujikawa, S.  Quintessence: A Review.  Class. Quant. Grav. 2013, 30, 214003, [arXiv:gr-qc/1304.1961].  doi:black10.1088/0264-9381/30/21/214003.
  * Hu et al. (2000) Hu, W.; Barkana, R.; Gruzinov, A.  Cold and fuzzy dark matter.  Phys. Rev. Lett. 2000, 85, 1158–1161, [astro-ph/0003365].  doi:black10.1103/PhysRevLett.85.1158.
  * Anagnostopoulos et al. (2019) Anagnostopoulos, F.K.; Benisty, D.; Basilakos, S.; Guendelman, E.I.  Dark energy and dark matter unification from dynamical space time: observational constraints and cosmological implications.  JCAP 2019, 1906, 003, [arXiv:gr-qc/1904.05762].  doi:black10.1088/1475-7516/2019/06/003.
  * Benisty et al. (2019) Benisty, D.; Guendelman, E.; Haba, Z.  Unification of dark energy and dark matter from diffusive cosmology.  Phys. Rev. 2019, D99, 123521, [arXiv:gr-qc/1812.06151].  doi:black10.1103/PhysRevD.99.123521.
  * Benisty and Guendelman (2018) Benisty, D.; Guendelman, E.I.  Unified dark energy and dark matter from dynamical spacetime.  Phys. Rev. 2018, D98, 023506, [arXiv:gr-qc/1802.07981].  doi:black10.1103/PhysRevD.98.023506.
  * Benisty and Guendelman (2017) Benisty, D.; Guendelman, E.I.  Interacting Diffusive Unified Dark Energy and Dark Matter from Scalar Fields.  Eur. Phys. J. 2017, C77, 396, [arXiv:gr-qc/1701.08667].  doi:black10.1140/epjc/s10052-017-4939-x.
  * Senovilla and Garfinkle (2015) Senovilla, J.M.M.; Garfinkle, D.  The 1965 Penrose singularity theorem.  Class. Quant. Grav. 2015, 32, 124008, [arXiv:gr-qc/1410.5226].  doi:black10.1088/0264-9381/32/12/124008.
  * De Felice and Tsujikawa (2010) De Felice, A.; Tsujikawa, S.  f(R) theories.  Living Rev. Rel. 2010, 13, 3, [arXiv:gr-qc/1002.4928].  doi:black10.12942/lrr-2010-3.
  * Capozziello and De Laurentis (2011) Capozziello, S.; De Laurentis, M.  Extended Theories of Gravity.  Phys. Rept. 2011, 509, 167–321, [arXiv:gr-qc/1108.6266].  doi:black10.1016/j.physrep.2011.09.003.
  * Nojiri et al. (2017) Nojiri, S.; Odintsov, S.D.; Oikonomou, V.K.  Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution.  Phys. Rept. 2017, 692, 1–104, [arXiv:gr-qc/1705.11098].  doi:black10.1016/j.physrep.2017.06.001.
  * Heisenberg (2019) Heisenberg, L.  A systematic approach to generalisations of General Relativity and their cosmological implications.  Phys. Rept. 2019, 796, 1–113, [arXiv:gr-qc/1807.01725].  doi:black10.1016/j.physrep.2018.11.006.
  * Bull et al. (2016) Bull, P.; others.  Beyond $\Lambda$CDM: Problems, solutions, and the road ahead.  Phys. Dark Univ. 2016, 12, 56–99, [arXiv:astro-ph.CO/1512.05356].  doi:black10.1016/j.dark.2016.02.001.
  * Banados and Ferreira (2010) Banados, M.; Ferreira, P.G.  Eddington’s theory of gravity and its progeny.  Phys. Rev. Lett. 2010, 105, 011101, [arXiv:astro-ph.CO/1006.1769].  [Erratum: Phys.Rev.Lett. 113, 119901 (2014)], doi:black10.1103/PhysRevLett.105.011101.
  * Alishahiha et al. (2004) Alishahiha, M.; Silverstein, E.; Tong, D.  DBI in the sky.  Phys. Rev. D 2004, 70, 123505, [hep-th/0404084].  doi:black10.1103/PhysRevD.70.123505.
  * Liu et al. (2012) Liu, Y.X.; Yang, K.; Guo, H.; Zhong, Y.  Domain Wall Brane in Eddington Inspired Born-Infeld Gravity.  Phys. Rev. D 2012, 85, 124053, [arXiv:hep-th/1203.2349].  doi:black10.1103/PhysRevD.85.124053.
  * Choudhury and Pal (2013) Choudhury, S.; Pal, S.  DBI Galileon inflation in background SUGRA.  Nucl. Phys. B 2013, 874, 85–114, [arXiv:hep-th/1208.4433].  doi:black10.1016/j.nuclphysb.2013.05.010.
  * Choudhury and Pal (2015) Choudhury, S.; Pal, S.  Primordial non-Gaussian features from DBI Galileon inflation.  Eur. Phys. J. C 2015, 75, 241, [arXiv:hep-th/1210.4478].  doi:black10.1140/epjc/s10052-015-3452-3.
  * Harko et al. (2013) Harko, T.; Lobo, F.S.N.; Mak, M.K.; Sushkov, S.V.  Structure of neutron, quark and exotic stars in Eddington-inspired Born-Infeld gravity.  Phys. Rev. D 2013, 88, 044032, [arXiv:gr-qc/1305.6770].  doi:black10.1103/PhysRevD.88.044032.
  * Wei et al. (2015) Wei, S.W.; Yang, K.; Liu, Y.X.  Black hole solution and strong gravitational lensing in Eddington-inspired Born–Infeld gravity.  Eur. Phys. J. C 2015, 75, 253, [arXiv:gr-qc/1405.2178].  [Erratum: Eur.Phys.J.C 75, 331 (2015)], doi:black10.1140/epjc/s10052-015-3556-9.
  * Shaikh (2015) Shaikh, R.  Lorentzian wormholes in Eddington-inspired Born-Infeld gravity.  Phys. Rev. D 2015, 92, 024015, [arXiv:gr-qc/1505.01314].  doi:black10.1103/PhysRevD.92.024015.
  * Avelino (2016) Avelino, P.P.  Inner Structure of Black Holes in Eddington-inspired Born-Infeld gravity: the role of mass inflation.  Phys. Rev. D 2016, 93, 044067, [arXiv:gr-qc/1511.03223].  doi:black10.1103/PhysRevD.93.044067.
  * Prasetyo et al. (2018) Prasetyo, I.; Husin, I.; Qauli, A.I.; Ramadhan, H.S.; Sulaksono, A.  Neutron stars in the braneworld within the Eddington-inspired Born-Infeld gravity.  JCAP 2018, 01, 027, [arXiv:astro-ph.CO/1708.04837].  doi:black10.1088/1475-7516/2018/01/027.
  * Chen et al. (2018) Chen, C.Y.; Bouhmadi-López, M.; Chen, P.  Black hole solutions in mimetic Born-Infeld gravity.  Eur. Phys. J. C 2018, 78, 59, [arXiv:gr-qc/1710.10638].  doi:black10.1140/epjc/s10052-018-5556-z.
  * Shaikh (2018) Shaikh, R.  Wormholes with nonexotic matter in Born-Infeld gravity.  Phys. Rev. D 2018, 98, 064033, [arXiv:gr-qc/1807.07941].  doi:black10.1103/PhysRevD.98.064033.
  * Jana et al. (2018) Jana, S.; Shaikh, R.; Sarkar, S.  Overcharging black holes and cosmic censorship in Born-Infeld gravity.  Phys. Rev. D 2018, 98, 124039, [arXiv:gr-qc/1808.09656].  doi:black10.1103/PhysRevD.98.124039.
  * Böhmer and Fiorini (2019) Böhmer, C.G.; Fiorini, F.  The regular black hole in four dimensional Born–Infeld gravity.  Class. Quant. Grav. 2019, 36, 12LT01, [arXiv:gr-qc/1901.02965].  doi:black10.1088/1361-6382/ab1e8d.
  * Delhom et al. (2019) Delhom, A.; Macedo, C.F.B.; Olmo, G.J.; Crispino, L.C.B.  Absorption by black hole remnants in metric-affine gravity.  Phys. Rev. 2019, D100, 024016, [arXiv:gr-qc/1906.06411].  doi:black10.1103/PhysRevD.100.024016.
  * Avelino and Ferreira (2012) Avelino, P.P.; Ferreira, R.Z.  Bouncing Eddington-inspired Born-Infeld cosmologies: an alternative to Inflation ?  Phys. Rev. D 2012, 86, 041501, [arXiv:astro-ph.CO/1205.6676].  doi:black10.1103/PhysRevD.86.041501.
  * Beltran Jimenez et al. (2018) Beltran Jimenez, J.; Heisenberg, L.; Olmo, G.J.; Rubiera-Garcia, D.  Born–Infeld inspired modifications of gravity.  Phys. Rept. 2018, 727, 1–129, [arXiv:gr-qc/1704.03351].  doi:black10.1016/j.physrep.2017.11.001.
  * Beltrán Jiménez and Delhom (2019) Beltrán Jiménez, J.; Delhom, A.  Ghosts in metric-affine higher order curvature gravity.  Eur. Phys. J. C 2019, 79, 656, [arXiv:gr-qc/1901.08988].  doi:black10.1140/epjc/s10052-019-7149-x.
  * Beltrán Jiménez and Delhom (2020) Beltrán Jiménez, J.; Delhom, A.  Instabilities in metric-affine theories of gravity with higher order curvature terms.  Eur. Phys. J. 2020, C80, 585, [arXiv:gr-qc/2004.11357].  doi:black10.1140/epjc/s10052-020-8143-z.
  * Delhom (2020) Delhom, A.  Minimal coupling in presence of non-metricity and torsion.  Eur. Phys. J. 2020, C80, 728, [arXiv:gr-qc/2002.02404].  doi:black10.1140/epjc/s10052-020-8330-y.
  * Afonso et al. (2018) Afonso, V.I.; Olmo, G.J.; Rubiera-Garcia, D.  Mapping Ricci-based theories of gravity into general relativity.  Phys. Rev. D 2018, 97, 021503, [arXiv:gr-qc/1801.10406].  doi:black10.1103/PhysRevD.97.021503.
  * Beltrán Jiménez et al. (2020) Beltrán Jiménez, J.; De Andrés, D.; Delhom, A.  Anisotropic deformations in a class of projectively-invariant metric-affine theories of gravity.  Class. Quant. Grav. 2020, 37, 225013, [arXiv:gr-qc/2006.07406].  doi:black10.1088/1361-6382/abb923.
  * Jana et al. (2018) Jana, S.; Chakravarty, G.K.; Mohanty, S.  Constraints on Born-Infeld gravity from the speed of gravitational waves after GW170817 and GRB 170817A.  Phys. Rev. D 2018, 97, 084011, [arXiv:gr-qc/1711.04137].  doi:black10.1103/PhysRevD.97.084011.
  * Barragan and Olmo (2010) Barragan, C.; Olmo, G.J.  Isotropic and Anisotropic Bouncing Cosmologies in Palatini Gravity.  Phys. Rev. 2010, D82, 084015, [arXiv:gr-qc/1005.4136].  doi:black10.1103/PhysRevD.82.084015.
  * Afonso et al. (2017) Afonso, V.I.; Olmo, G.J.; Rubiera-Garcia, D.  Scalar geons in Born-Infeld gravity.  JCAP 2017, 08, 031, [arXiv:gr-qc/1705.01065].  doi:black10.1088/1475-7516/2017/08/031.
  * Solà et al. (2017) Solà, J.; Gómez-Valent, A.; de Cruz Pérez, J.  First evidence of running cosmic vacuum: challenging the concordance model.  Astrophys. J. 2017, 836, 43, [arXiv:astro-ph.CO/1602.02103].  doi:black10.3847/1538-4357/836/1/43.
  * Barrow et al. (2021) Barrow, J.D.; Basilakos, S.; Saridakis, E.N.  Big Bang Nucleosynthesis constraints on Barrow entropy.  Phys. Lett. 2021, B815, 136134, [arXiv:gr-qc/2010.00986].  doi:black10.1016/j.physletb.2021.136134.
  * Jiménez et al. (2021) Jiménez, J.B.; Delhom, A.; Olmo, G.J.; Orazi, E.  Born-Infeld gravity: Constraints from light-by-light scattering and an effective field theory perspective.  Phys. Lett. B 2021, 820, 136479, [arXiv:gr-qc/2104.01647].  doi:black10.1016/j.physletb.2021.136479.
  * Latorre et al. (2018) Latorre, A.D.I.; Olmo, G.J.; Ronco, M.  Observable traces of non-metricity: new constraints on metric-affine gravity.  Phys. Lett. 2018, B780, 294–299, [arXiv:hep-th/1709.04249].  doi:black10.1016/j.physletb.2018.03.002.
  * Delhom et al. (2020) Delhom, A.; Miralles, V.; Peñuelas, A.  Effective interactions in Ricci-Based Gravity below the non-metricity scale.  Eur. Phys. J. 2020, C80, 340, [arXiv:hep-th/1907.05615].  doi:black10.1140/epjc/s10052-020-7880-3.
  * Jimenez and Loeb (2002) Jimenez, R.; Loeb, A.  Constraining cosmological parameters based on relative galaxy ages.  Astrophys. J. 2002, 573, 37–42, [arXiv:astro-ph/astro-ph/0106145].  doi:black10.1086/340549.
  * Moresco et al. (2012a) Moresco, M.; Verde, L.; Pozzetti, L.; Jimenez, R.; Cimatti, A.  New constraints on cosmological parameters and neutrino properties using the expansion rate of the Universe to z 1.75.  JCAP 2012, 1207, 053, [arXiv:astro-ph.CO/1201.6658].  doi:black10.1088/1475-7516/2012/07/053.
  * Moresco et al. (2012b) Moresco, M.; others.  Improved constraints on the expansion rate of the Universe up to z 1.1 from the spectroscopic evolution of cosmic chronometers.  JCAP 2012, 1208, 006, [arXiv:astro-ph.CO/1201.3609].  doi:black10.1088/1475-7516/2012/08/006.
  * Moresco (2015) Moresco, M.  Raising the bar: new constraints on the Hubble parameter with cosmic chronometers at z 2.  Mon. Not. Roy. Astron. Soc. 2015, 450, L16–L20, [arXiv:astro-ph.CO/1503.01116].  doi:black10.1093/mnrasl/slv037.
  * Moresco et al. (2016) Moresco, M.; Pozzetti, L.; Cimatti, A.; Jimenez, R.; Maraston, C.; Verde, L.; Thomas, D.; Citro, A.; Tojeiro, R.; Wilkinson, D.  A 6% measurement of the Hubble parameter at $z\sim 0.45$: direct evidence of the epoch of cosmic re-acceleration.  JCAP 2016, 1605, 014, [arXiv:astro-ph.CO/1601.01701].  doi:black10.1088/1475-7516/2016/05/014.
  * Scolnic et al. (2018) Scolnic, D.; others.  The Complete Light-curve Sample of Spectroscopically Confirmed SNe Ia from Pan-STARRS1 and Cosmological Constraints from the Combined Pantheon Sample.  Astrophys. J. 2018, 859, 101, [arXiv:astro-ph.CO/1710.00845].  doi:black10.3847/1538-4357/aab9bb.
  * Anagnostopoulos et al. (2020) Anagnostopoulos, F.K.; Basilakos, S.; Saridakis, E.N.  Observational constraints on Barrow holographic dark energy.  Eur. Phys. J. 2020, C80, 826, [arXiv:gr-qc/2005.10302].  doi:black10.1140/epjc/s10052-020-8360-5.
  * Roberts et al. (2017) Roberts, C.; Horne, K.; Hodson, A.O.; Leggat, A.D.  Tests of $\Lambda$CDM and Conformal Gravity using GRB and Quasars as Standard Candles out to $z\sim 8$ 2017.  [arXiv:astro-ph.CO/1711.10369].
  * Demianski et al. (2017) Demianski, M.; Piedipalumbo, E.; Sawant, D.; Amati, L.  Cosmology with gamma-ray bursts: I. The Hubble diagram through the calibrated $E_{\rm p,i}$ \- $E_{\rm iso}$ correlation.  Astron. Astrophys. 2017, 598, A112, [arXiv:astro-ph.CO/1610.00854].  doi:black10.1051/0004-6361/201628909.
  * Hogg et al. (2020) Hogg, N.B.; Martinelli, M.; Nesseris, S.  Constraints on the distance duality relation with standard sirens 2020.  [arXiv:astro-ph.CO/2007.14335].
  * Martinelli et al. (2020) Martinelli, M.; others.  Euclid: Forecast constraints on the cosmic distance duality relation with complementary external probes 2020.  [arXiv:astro-ph.CO/2007.16153].
  * Benisty and Staicova (2020) Benisty, D.; Staicova, D.  Testing Low-Redshift Cosmic Acceleration with the Complete Baryon Acoustic Oscillations data collection 2020.  [arXiv:astro-ph.CO/2009.10701].
  * Percival et al. (2010) Percival, W.J.; others.  Baryon Acoustic Oscillations in the Sloan Digital Sky Survey Data Release 7 Galaxy Sample.  Mon. Not. Roy. Astron. Soc. 2010, 401, 2148–2168, [arXiv:astro-ph.CO/0907.1660].  doi:black10.1111/j.1365-2966.2009.15812.x.
  * Beutler et al. (2011) Beutler, F.; Blake, C.; Colless, M.; Jones, D.H.; Staveley-Smith, L.; Campbell, L.; Parker, Q.; Saunders, W.; Watson, F.  The 6dF Galaxy Survey: Baryon Acoustic Oscillations and the Local Hubble Constant.  Mon. Not. Roy. Astron. Soc. 2011, 416, 3017–3032, [arXiv:astro-ph.CO/1106.3366].  doi:black10.1111/j.1365-2966.2011.19250.x.
  * Busca et al. (2013) Busca, N.G.; others.  Baryon Acoustic Oscillations in the Ly-$\alpha$ forest of BOSS quasars.  Astron. Astrophys. 2013, 552, A96, [arXiv:astro-ph.CO/1211.2616].  doi:black10.1051/0004-6361/201220724.
  * Anderson et al. (2013) Anderson, L.; others.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations in the Data Release 9 Spectroscopic Galaxy Sample.  Mon. Not. Roy. Astron. Soc. 2013, 427, 3435–3467, [arXiv:astro-ph.CO/1203.6594].  doi:black10.1111/j.1365-2966.2012.22066.x.
  * Seo et al. (2012) Seo, H.J.; others.  Acoustic scale from the angular power spectra of SDSS-III DR8 photometric luminous galaxies.  Astrophys. J. 2012, 761, 13, [arXiv:astro-ph.CO/1201.2172].  doi:black10.1088/0004-637X/761/1/13.
  * Ross et al. (2015) Ross, A.J.; Samushia, L.; Howlett, C.; Percival, W.J.; Burden, A.; Manera, M.  The clustering of the SDSS DR7 main Galaxy sample – I. A 4 per cent distance measure at $z=0.15$.  Mon. Not. Roy. Astron. Soc. 2015, 449, 835–847, [arXiv:astro-ph.CO/1409.3242].  doi:black10.1093/mnras/stv154.
  * Tojeiro et al. (2014) Tojeiro, R.; others.  The clustering of galaxies in the SDSS-III Baryon Oscillation Spectroscopic Survey: galaxy clustering measurements in the low redshift sample of Data Release 11.  Mon. Not. Roy. Astron. Soc. 2014, 440, 2222–2237, [arXiv:astro-ph.CO/1401.1768].  doi:black10.1093/mnras/stu371.
  * Bautista et al. (2018) Bautista, J.E.; others.  The SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations at redshift of 0.72 with the DR14 Luminous Red Galaxy Sample.  Astrophys. J. 2018, 863, 110, [arXiv:astro-ph.CO/1712.08064].  doi:black10.3847/1538-4357/aacea5.
  * de Carvalho et al. (2018) de Carvalho, E.; Bernui, A.; Carvalho, G.C.; Novaes, C.P.; Xavier, H.S.  Angular Baryon Acoustic Oscillation measure at $z=2.225$ from the SDSS quasar survey.  JCAP 2018, 1804, 064, [arXiv:astro-ph.CO/1709.00113].  doi:black10.1088/1475-7516/2018/04/064.
  * Ata et al. (2018) Ata, M.; others.  The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey DR14 quasar sample: first measurement of baryon acoustic oscillations between redshift 0.8 and 2.2.  Mon. Not. Roy. Astron. Soc. 2018, 473, 4773–4794, [arXiv:astro-ph.CO/1705.06373].  doi:black10.1093/mnras/stx2630.
  * Abbott et al. (2019) Abbott, T.M.C.; others.  Dark Energy Survey Year 1 Results: Measurement of the Baryon Acoustic Oscillation scale in the distribution of galaxies to redshift 1.  Mon. Not. Roy. Astron. Soc. 2019, 483, 4866–4883, [arXiv:astro-ph.CO/1712.06209].  doi:black10.1093/mnras/sty3351.
  * Molavi and Khodam-Mohammadi (2019) Molavi, Z.; Khodam-Mohammadi, A.  Observational tests of Gauss-Bonnet like dark energy model.  Eur. Phys. J. Plus 2019, 134, 254, [arXiv:gr-qc/1906.05668].  doi:black10.1140/epjp/i2019-12723-x.
  * Riess et al. (2021) Riess, A.G.; Casertano, S.; Yuan, W.; Bowers, J.B.; Macri, L.; Zinn, J.C.; Scolnic, D.  Cosmic Distances Calibrated to 1% Precision with Gaia EDR3 Parallaxes and Hubble Space Telescope Photometry of 75 Milky Way Cepheids Confirm Tension with $\Lambda$CDM.  Astrophys. J. Lett. 2021, 908, L6, [arXiv:astro-ph.CO/2012.08534].  doi:black10.3847/2041-8213/abdbaf.
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D.; Hogg, D.W.; Lang, D.; Goodman, J.  emcee: The MCMC Hammer.  Publ. Astron. Soc. Pac. 2013, 125, 306–312, [arXiv:astro-ph.IM/1202.3665].  doi:black10.1086/670067.
  * Handley et al. (2015) Handley, W.J.; Hobson, M.P.; Lasenby, A.N.  PolyChord: nested sampling for cosmology.  Mon. Not. Roy. Astron. Soc. 2015, 450, L61–L65, [arXiv:astro-ph.CO/1502.01856].  doi:black10.1093/mnrasl/slv047.
  * Lewis (2019) Lewis, A.  GetDist: a Python package for analysing Monte Carlo samples 2019.  [arXiv:astro-ph.IM/1910.13970].
  * Beltran Jimenez et al. (2017) Beltran Jimenez, J.; Heisenberg, L.; Olmo, G.J.; Rubiera-Garcia, D.  On gravitational waves in Born-Infeld inspired non-singular cosmologies.  JCAP 2017, 10, 029, [arXiv:hep-th/1707.08953].  [Erratum: JCAP 08, E01 (2018)], doi:black10.1088/1475-7516/2017/10/029.
  * Escamilla-Rivera et al. (2012) Escamilla-Rivera, C.; Banados, M.; Ferreira, P.G.  A tensor instability in the Eddington inspired Born-Infeld Theory of Gravity.  Phys. Rev. D 2012, 85, 087302, [arXiv:gr-qc/1204.1691].  doi:black10.1103/PhysRevD.85.087302.