# Covariant formulation of refracted gravity

Andrea P. Sannaa,b,c<EMAIL_ADDRESS>Titos Matsakosc
<EMAIL_ADDRESS>Antonaldo Diaferioc,d<EMAIL_ADDRESS>aDipartimento di Fisica, Università di Cagliari, Cittadella Universitaria,
09042, Monserrato, Italy  bIstituto Nazionale di Fisica Nucleare (INFN),
Sezione di Cagliari, Cittadella Universitaria, 09042, Monserrato, Italy
cDipartimento di Fisica, Università di Torino, via P. Giuria 1, 10125, Torino,
Italy  dIstituto Nazionale di Fisica Nucleare (INFN), Sezione di Torino,
Torino, Italy

(August 27, 2024)

###### Abstract

We propose a covariant formulation of refracted gravity (RG), a classical
theory of gravity based on the introduction of the gravitational permittivity
— a monotonic function of the local mass density — in the standard Poisson
equation. The gravitational permittivity mimics the dark matter phenomenology.
Our covariant formulation of RG (CRG) belongs to the class of scalar-tensor
theories, where the scalar field $\varphi$ has a self-interaction potential
$\mathcal{V}(\varphi)=-\Xi\varphi$, with $\Xi$ a normalization constant. We
show that the scalar field is twice the gravitational permittivity in the
weak-field limit. Far from a spherical source of density $\rho_{s}(r)$, the
transition between the Newtonian and the RG regime appears below the
acceleration scale $a_{\Xi}=(2\Xi-8\pi G\rho/\varphi)^{1/2}$, with
$\rho=\rho_{s}+\rho_{\mathrm{bg}}$ and $\rho_{\mathrm{bg}}$ an isotropic and
homogeneous background. In the limit $2\Xi\gg 8\pi G\rho/\varphi$, we obtain
$a_{\Xi}\sim 10^{-10}$ m s-2. This acceleration is comparable to the
acceleration $a_{0}$ originally introduced in Modified Newtonian Dynamics
(MOND). From CRG, we also derive the modified Friedmann equations for an
expanding, homogeneous, and isotropic universe. We find that the same scalar
field that mimics dark matter also drives the accelerated expansion of the
Universe. Since $\Xi$ plays a role roughly similar to the cosmological
constant $\Lambda$ in the standard model and has a comparable value, CRG
suggests a natural explanation of the known relation $a_{0}\sim\Lambda^{1/2}$.
CRG thus appears to describe both the dynamics of cosmic structure and the
expanding Universe with a single scalar field, and falls within the family of
models that unify the two dark sectors, highlighting a possible deep
connection between phenomena currently attributed to dark matter and dark
energy separately.

##  I Introduction

The current cosmological model, $\Lambda$CDM, assumes that gravitational
interactions are ruled by General Relativity (GR); the model relies on the
existence of non-baryonic Cold Dark Matter (CDM), and a positive cosmological
constant $\Lambda$ 1995Natur.377..600O . Non-baryonic dark matter is required
to account for the abundance of light elements 2016RvMP…88a5004C , the
amplitude of the power spectrum of the temperature anisotropies in the Cosmic
Microwave Background (CMB) Aghanim:2019ame , and the dynamics of cosmic
structures Clowe:2006eq ; Dodelson:2001ux ; DelPopolo:2013qba ; Akrami:2018mcd
; Planck:2018nkj . The cosmological density parameter $\Omega_{\Lambda
0}=\Lambda/3H_{0}^{2}$ associated to $\Lambda$, with $H_{0}$ the Hubble
constant, accounts for the negative value of the deceleration parameter
$q_{0}$ measured with the Hubble diagram of type Ia supernovae Riess+98 ;
Perlmutter+99 . The curvature of the Universe, $\Omega_{k}=\Omega_{\Lambda
0}+\Omega_{0}-1$, measured from the CMB power spectrum, suggests a null
curvature Aghanim:2018eyx .

Although the $\Lambda$CDM model agrees with most of the rich observational
information currently available, a number of tensions are present both on
large cosmic scales and on the scale of galaxies. The Hubble constant $H_{0}$
estimated with the distance ladder in the local Universe 2016ApJ…826…56R ;
2019ApJ…876…85R  is more than $4\sigma$ larger than $H_{0}$ inferred from the
measurements of the CMB Verde:2019ivm . The observed abundances of the light
elements are consistent with the Big Bang nucleosynthesis scenario, except
7Li, whose abundance is a factor of $\sim 3$ smaller Mathews:2019hbi .
Numerous features of the CMB temperature anistropies are present on large
scales Ade:2015hxq . The probability of some of these features to appear
individually is $\lesssim 0.1\%$; the combined probability of the uncorrelated
features is $\lesssim 0.001\%$ and might represent a serious challenge to the
$\Lambda$CDM model Schwarz:2015cma . The lensing amplitude in the CMB power
spectra is enhanced compared to $\Lambda$CDM expectations Aghanim:2019ame  and
would suggest a positive rather than a null curvature of the Universe
DiValentino:2019qzk . A slight tension also appears for the normalization of
the power spectrum $\sigma_{8}$ (e.g. Raveri:2015maa ): the value inferred
from the CMB measurements Aghanim:2018eyx  is $2\sigma$ larger than the value
derived from the tomographic weak gravitational lensing analysis of the KiDS
imaging data Hildebrandt:2016iqg .

In addition, the cosmological constant poses a fine-tuning problem that is
theoretically challenging. If we associate $\Lambda$ with the ground state
energy level of the vacuum in quantum field theory, its measured value,
$\Lambda\sim 10^{-12}$ eV4, appears to be $\sim 120$ orders of magnitude
smaller than expected Weinberg+89 ; Padilla:2015aaa . If we associate
$\Lambda$ with the energy scale up to which the Standard Model of particle
physics has been tested, $\sim 1$ TeV, the discrepancy reduces to $\sim 60$
orders of magnitude Joyce+15 , but it remains severe. The most popular
solution to the $\Lambda$ problem is to suppress $\Lambda$ in the Einstein-
Hilbert action and generate the accelerated expansion of the Universe with
dark energy, an auxiliary scalar field with proper features. The specific
implementation of this idea has generated a large number of different models,
that may or may not modify Einstein’s equations (see Peebles+03 ;
Copeland:2006wr ; Joyce+15 ; Amendola:2016saw  for extensive reviews).

Additionally, dark matter poses an astrophysical and an astro-particle
problem. On the scale of galaxies, a number of issues, including the core-cusp
problem in dwarf galaxies, the missing satellite problem, and the plane of
satellite galaxies in the Milky Way and nearby large galaxies, are hardly
explained within the CDM framework (e.g. DelPopolo+17 ; Kroupa2012 ), and
some, like the radial acceleration relation or the baryonic Tully-Fisher
relation in disk galaxies, would require interactions between dark and
baryonic matter that should be unrealistically fine-tuned (see, e.g.
deMartino:2020gfi  for a review).

Weakly interacting massive particles (WIMPs) are the most plausible candidate
to make up the non-baryonic dark matter Bertone+05 . However, attempts to,
directly or indirectly, detect these particles have not yet produced
unquestionable results Tanabashi:2018oca . A possible solution to this impasse
is to suppose that the dark matter is made up of more exotic particles, for
example sterile neutrinos, axions, or ultra-light bosons. The detection of
these particles requires direct or indirect experiments different from those
conceived for detecting WIMPs Buckley:2017ijx . These exotic particles might
also solve some, albeit not all, of the problems on the scale of galaxies
(e.g. deMartino:2020gfi ). The phenomenology of galaxies appears to be
reproduced if the dark matter particles have distinct peculiar features, like
superfluidity Berezhiani:2015bqa  or gravitational polarization
Blanchet:2006yt ; Blanchet:2017duj .

Alternatively, the dynamics on the scale of galaxies could be explained by
modifying the theory of gravity in the weak-field Newtonian limit without
resorting to the existence of dark matter. The idea of a Modified Newtonian
Dynamics (MOND) was originally motivated only by the observations of the flat
rotation curves of disk galaxies Milgrom1 , and many of the observational
challenges on the scale of galaxies were actually predicted by MOND Milgrom2 ;
Milgrom3 ; Sanders:2002pf . The MOND formulation is purely phenomenological
and its extension to a covariant formulation has proven to be challenging
Famaey+12 ; MilgromReview ; Skordis:2019fxt ; Zlosnik:2017xpr .

Although the problems of the cosmological constant and dark matter are usually
considered two separate issues, attempts to unify the two dark sectors in a
single framework are numerous. For example, Ferreira et al. Ferreira:2018wup
suggest a model where the dark matter is made of two superfluids arising from
two distinguishable states separated by a small energy difference: particles
of one species can be converted into the other and the interaction between
these multi-state components of the dark matter can drive cosmic acceleration.
Alternatively, in emergent gravity, where both the classical space-time
structure and gravity emerge from an underlying microscopic quantum theory
Sakharov:1967pk ; Padmanabhan:2014jta ; Verlinde:2016toy , the two dark
sectors can be unified when, for example, the dark energy fluid is modelled as
a critical Bose-Einstein condensate of gravitons Cadoni:2017evg ;
Cadoni:2018dnd ; Tuveri:2019zor ; Cadoni:2020izk .

Attempts to describe both the accelerated expansion of the Universe and the
dynamics of cosmic structures without dark matter by building a modified
theory of gravity are also numerous (see e.g. Clifton+12 ). According to
Lovelock’s theorem Lovelock:1971yv ; Lovelock:1972vz , we can build a metric
theory of gravity different from GR by, for example, allowing derivatives of
the metric tensor higher than second order in the field equations, or
introducing other fields in addition to the metric tensor.

Conformal gravity adopts the former route and replaces the Einstein-Hilbert
action with the contraction of the fourth-rank conformal tensor introduced by
Weyl Mannheim:1992tr . Conformal gravity does not present ghost-like
instabilities, that might be common in theories with high-order derivatives
Bender:2007wu , and is renormalizable Mannheim:2011ds . Unfortunately,
although conformal gravity successfully reproduces the accelerated expansion
of the Universe Diaferio:2011kc , the expected abundance of primordial
deuterium is orders of magnitudes smaller than observed Elizondo:1994vh .
Moreover, conformal gravity is unable to reproduce the dynamics of galaxy
clusters 8219260 ; Diaferio:2008gh , and appears to require a fine-tuning
condition to describe the phenomenology of gravitational lensing and the
dynamics of disk galaxies 2019CQGra..36x5014C .

We can preserve the second-order field equations by introducing a single
scalar field that drives both the accelerated expansion of the Universe and
the formation of cosmic structure (e.g. Carneiro:2018zet ). The case of a
classical scalar field with a non-canonical kinetic term in its Lagrangian
generates the class of Unified Dark Matter models Bertacca:2010ct . These
models can be a viable alternative to $\Lambda$CDM Camera:2010wm ;
Camera:2017tws  if the effective sound speed is small enough that the scalar
field can cluster Bertacca:2008uf ; Camera:2009uz .

Here, we contribute to the quest for a unified model of the two dark sectors
by proposing a novel scalar-tensor theory Fujiibook ; Quiros2019  where the
scalar field is responsible for both the dynamics of cosmic structures and the
accelerated expansion of the Universe. This scalar-tensor theory is the
covariant formulation of refracted gravity (RG), a new phenomenological
modified theory of gravity proposed by Matsakos & Diaferio MatsakosDiaferio16
. RG appears to reproduce the phenomenology on the scale of galaxies and
galaxy clusters by introducing a monotonic function of the local mass-density
in the standard Poisson equation, termed gravitational permittivity. Indeed,
Cesare et al. Cesare+20  showed that RG properly describes the rotation curves
and the vertical velocity dispersion profiles of 30 disk galaxies in the
DiskMass Survey Bershady:2010sj , and the dynamics of stars and globular
clusters in the outer regions of three elliptical galaxies of type E0
Cesare+21 . Here, we provide a covariant formulation of this non-relativistic
formulation of RG.

Sect. II reviews the relevant features of RG. In Sect. III, we derive the
covariant formulation of RG (CRG, hereafter) in the framework of scalar-tensor
theories, and show that the scalar field can be identified with the
gravitational permittivity. In Sect. IV, we consider a homogeneous and
isotropic universe and derive the modified Friedmann equations; we show that
the scalar field is responsible for the accelerated expansion of the Universe
and that CRG can properly describe the Hubble diagram of the high-redsfhit
type Ia supernovae. We conclude in Sect. V.

##  II Non-relativistic Refracted Gravity

The phenomenological RG is based on the modified Poisson equation (hereafter
RG equation) MatsakosDiaferio16 ,

$\vec{\nabla}\cdot\left(\epsilon\vec{\nabla}\Phi\right)=4\pi G\rho\,,$ (1)

where $\Phi$ is the gravitational potential, $G$ the gravitational constant,
and $\rho$ the density of ordinary matter. The parameter $\epsilon$ is the
gravitational permittivity, according to the mathematical similarity with the
term on the left-hand side of the Poisson equation that describes electric
fields in matter. We emphasize that no other parallels have been drawn with
electrodynamics. As a starting hypothesis, the permittivity was prescribed to
be a monotonic function of the local mass density, i.e. $\epsilon(\rho)$.
However, $\rho$ was only chosen because is the simplest scalar field
characterizing the matter distribution in the weak-field regime; $\epsilon$
could in principle depend on other local quantities, for example the trace of
the energy-momentum tensor or the entropy.

On scales where the visible mass can accurately explain the observed dynamics
according to Newtonian gravity, for example on the scale of stars, $\epsilon$
must be constant and equal to $1$ in order to recover the standard Newtonian
Poisson equation, $\nabla^{2}\Phi=4\pi G\rho$. Therefore, we can adopt the
form of the permittivity

$\epsilon(\rho)\simeq\left\\{\begin{array}[]{ll}1&\mathrm{for}\ \
\rho\gtrsim\rho_{\mathrm{thr}}\\\ \epsilon_{\mathrm{v}}&\mathrm{for}\ \
\rho\ll\rho_{\mathrm{thr}}\end{array}\right.\,,$ (2)

where $\rho_{\mathrm{thr}}$ is the density threshold that sets the transition
between the Newtonian regime and the modified gravity regime;111Matsakos &
Diaferio MatsakosDiaferio16  adopt the designation critical density
$\rho_{\mathrm{c}}$ rather than threshold density $\rho_{\mathrm{thr}}$. Here,
we prefer to adopt the latter, to avoid any confusion with the critical
density of the Universe, used in modelling a homogeneous and isotropic
universe. In addition, here $\rho_{\mathrm{thr}}$ is explicitly related to the
value below which $\epsilon$ deviates from unity. $0<\epsilon_{\mathrm{v}}<1$
is the permittivity of the vacuum.

The presence of $\epsilon$ in the Poisson equation, Eq. (1), has two main
effects: (1) in low-density regions, it generates a stronger gravitational
field than in Newtonian gravity; (2) it bends the gravitational field lines in
regions where $\vec{\nabla}\Phi$ and $\vec{\nabla}\epsilon$ are not parallel.

The former effect is trivially seen in the simple case
$\epsilon=\mathrm{const}<1$, when $\rho<\rho_{\mathrm{thr}}$: Eq. (1) becomes
$\nabla^{2}\Phi=4\pi G\rho/\epsilon$, and shows that the resulting field is
equivalent to a Newtonian field originating from a larger “effective” mass
density $\rho/\epsilon$, or a larger gravitational constant $G/\epsilon$. When
$\vec{\nabla}\Phi$ and $\vec{\nabla}\epsilon$ are not parallel, the field
lines are refracted; in other words, the field lines change their direction
when they cross the iso-surfaces of $\epsilon$ at an angle different from
$\pi/2$. This effect can also generate non-Newtonian phenomena in regions
where the density is larger than $\rho_{\mathrm{thr}}$, because the
redirection of the field lines has non-local consequences. For example, when
we consider only ordinary matter as the gravitational source, RG predicts flat
rotation curves in disk galaxies even in regions where
$\rho>\rho_{\mathrm{thr}}$ Cesare+20 . The redirection of the field lines is
also expected to be responsible for the mass discrepancy in galaxy clusters
MatsakosDiaferio16 . In these two studies, $\epsilon(\rho)$ is assumed to be a
smooth step function, with the density threshold $\rho_{\mathrm{thr}}$ in the
range $10^{-27}$–$10^{-24}$ g cm-3.

##  III Refracted Gravity as a scalar-tensor theory

The family of scalar-tensor (TeS) theories, with the scalar field, $\varphi$,
non-minimally coupled to the rank-2 tensor field $g_{\mu\nu}$, can be derived
from the action222Hereafter, we use natural units with $c=\hbar=1$.

$\displaystyle\mathcal{S}$ $\displaystyle=\frac{1}{16\pi G}\int d^{4}x\
\sqrt{g}\left[\varphi
R+\frac{\mathcal{W}(\varphi)}{\varphi}\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi+2\mathcal{V}(\varphi)\right]$
$\displaystyle\quad+\int d^{4}x\sqrt{g}\
\mathcal{L}_{\mathrm{m}}(g_{\mu\nu},\psi_{\mathrm{m}})\,,$ (3)

where $-g$ is the determinant of $g_{\mu\nu}$, $R=g^{\mu\nu}R_{\mu\nu}$ the
Ricci scalar, $R_{\mu\nu}$ the Ricci tensor, $g^{\mu\nu}$ the inverse metric,
$\nabla_{\mu}$ the covariant derivative,
$\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi\equiv
g^{\alpha\beta}\nabla_{\alpha}\varphi\nabla_{\beta}\varphi$, and
$\mathcal{L}_{\mathrm{m}}$ the matter Lagrangian density, with
$\psi_{\mathrm{m}}$ the matter fields Faraoni+04 . The potential,
$\mathcal{V}(\varphi)$, and the general differentiable function of the scalar
field, $\mathcal{W}(\varphi)$, parametrize the family of TeS theories. When
$\varphi\to 1$, the TeS action reduces to the standard Einstein-Hilbert action
with a cosmological constant equal to $-\mathcal{V}(1)$. Hereafter, for the
sake of simplicity, we will consider the $\varphi$-dependence of $\mathcal{W}$
and $\mathcal{V}$ implicit.

Varying Eq. (3) with respect to the metric yields the ten modified Einstein
field equations in the Jordan frame (JF)333We adopt the sign conventions of
Weinberg weinberggravitation , with the Ricci tensor given by
$R_{\mu\nu}=\partial_{\nu}\Gamma^{\alpha}_{\mu\alpha}-\partial_{\alpha}\Gamma^{\alpha}_{\mu\nu}+\Gamma^{\alpha}_{\mu\beta}\Gamma^{\beta}_{\alpha\nu}-\Gamma^{\alpha}_{\mu\nu}\Gamma^{\beta}_{\alpha\beta}$,
and the standard Einstein equations $G_{\mu\nu}=-8\pi GT_{\mu\nu}$.

$\displaystyle\varphi\left(R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\right)-\left(\frac{\mathcal{W}}{2\varphi}\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi+\Box\varphi+\mathcal{V}\right)g_{\mu\nu}$
$\displaystyle\quad+\nabla_{\mu}\nabla_{\nu}\varphi+\frac{\mathcal{W}}{\varphi}\nabla_{\mu}\varphi\nabla_{\nu}\varphi=-8\pi
GT_{\mu\nu}\,,$ (4)

where $\Box\equiv g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the d’Alembertian
operator (e.g., Quiros+16 ). Varying Eq. (3) with respect to the scalar field
yields the equation for $\varphi$

$\varphi
R+\left(\frac{\mathcal{W}}{\varphi}-\frac{d\mathcal{W}}{d\varphi}\right)\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi-2\mathcal{W}\Box\varphi+2\varphi\frac{d\mathcal{V}}{d\varphi}=0\,.$
(5)

We define CRG by setting

$\mathcal{W}(\varphi)=-1\,.$ (6)

We note that the original Brans-Dicke theory — with a constant $\mathcal{W}$
of order unity and a zero potential — is ruled out by solar system
experimental tests, which give the constraint $\mathcal{W}\gtrsim 40\,000$
Faraoni+04 ; Clifton+12 , and by recent results from CMB observations (e.g.,
Li2013 ; Avilez2014 ). However, these constraints can be avoided by adding a
self-interaction potential Hrycyna2014 ; Quiros2019 , which we define as

$\mathcal{V}(\varphi)=-\Xi\varphi\,,$ (7)

with $\Xi$ a constant parameter.

With the definitions of Eqs. (6) and (7), the modified field equations and the
equation for $\varphi$, in the JF, become

$\displaystyle\varphi\left(R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\right)-\left(-\frac{1}{2\varphi}\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi+\Box\varphi-\Xi\varphi\right)g_{\mu\nu}$
$\displaystyle\quad+\nabla_{\mu}\nabla_{\nu}\varphi-\frac{1}{\varphi}\nabla_{\mu}\varphi\nabla_{\nu}\varphi=-8\pi
GT_{\mu\nu}\,,$ (8) $\displaystyle\varphi
R-\frac{1}{\varphi}\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi+2\Box\varphi-2\Xi\varphi=0\,.$
(9)

By using Eq. (9) to simplify Eq. (8), and by using Eq. (8), contracted with
$g^{\mu\nu}$, to simplify the resulting Eq. (9), we obtain the following CRG
equations:

$\displaystyle\varphi
R_{\mu\nu}+\nabla_{\mu}\nabla_{\nu}\varphi-\frac{1}{\varphi}\nabla_{\mu}\varphi\nabla_{\nu}\varphi$
$\displaystyle=-8\pi GT_{\mu\nu}\,,$ (10)
$\displaystyle\Box\varphi-2\Xi\varphi$ $\displaystyle=8\pi GT\,.$ (11)

###  III.1 Relation to Horndeski theories

TeS theories belong to the wider family of Horndeski models (for recent
reviews on the topic, see, e.g., Joyce+15 ; Kobayashi+19 ), and, consequently,
CRG can also be identified as a special case of Horndeski theories. The
detection of GW170817 and GRB170817A Kase:2018aps ; Wang:2017rpx ;
Sakstein:2017xjx  has greatly constrained the Horndeski theory space. The
currently allowed Lagrangian is (Noller+19 ; Creminelli+17 ; Sakstein+17  and
references therein)

$\mathcal{L}_{\mathrm{H}}=\mathcal{G}_{4}\left(\varphi\right)R+\mathcal{G}_{2}\left(\varphi,X\right)-\mathcal{G}_{3}\left(\varphi,X\right)\Box\varphi\,,$
(12)

where $\mathcal{G}_{i}$’s are functions of the scalar field $\varphi$ and
$X\equiv\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi$ is the kinetic term. By
comparing Eq. (12) with the Lagrangian of Eq. (3), and by ignoring the factor
$1/16\pi G$, we infer that $\mathcal{G}_{4}\leftrightarrow\varphi$,
$\mathcal{G}_{2}\leftrightarrow-\frac{1}{\varphi}\nabla^{\alpha}\varphi\nabla_{\alpha}\varphi-2\Xi\varphi$,
and $\mathcal{G}_{3}\leftrightarrow 0$.444The $f(R)$ models of gravity are
also obtained from the Lagrangian of Eq. (12) by setting $\mathcal{G}_{3}=0$.
This result is expected because of the equivalence between the Brans-Dicke
theory and $f(R)$ gravity f(R)Faraoni . Therefore, CRG is expected to inherit
the important properties of the viable Horndeski theories, including the
absence of ghost-like degrees of freedom, i.e. the Ostrogradski instabilities
Gleyzes+14 .

###  III.2 Screening mechanism of CRG

GR accurately describes the gravitational interactions on the scale of stars
and smaller scales, including the strong gravitational field regime.
Therefore, any attempt to modify the theory of gravity by adding new degrees
of freedom must provide a screening mechanism to avoid detectable
discrepancies in the local tests of gravity Joyce+15 .

The screening mechanism depends on the local mass density and/or the local
gravitational potential. It is thus convenient to study the screening
mechanism in the Einstein frame (EF), where the scalar field is minimally
coupled to both gravity and the matter fields. The transition from the JF to
the EF can be performed with the conformal transformation of the metric
Clifton+12

$\tilde{g}_{\mu\nu}=e^{-2\Gamma}g_{\mu\nu}\quad\mathrm{and}\quad\sqrt{\tilde{g}}=e^{-4\Gamma}\sqrt{g}\,,$
(13)

where the tilde indicates quantities computed in the EF and $\Gamma$ is a
function of the scalar field.

From the JF action of Eq. (3), the corresponding action in the EF is

$\displaystyle\tilde{\mathcal{S}}$ $\displaystyle=\frac{1}{4\pi G}\int d^{4}x\
\sqrt{\tilde{g}}\left[\frac{\tilde{R}}{4}-\frac{1}{2}\tilde{g}^{\mu\nu}\tilde{\nabla}_{\mu}\Psi\tilde{\nabla}_{\nu}\Psi-\mathcal{U}(\Psi)\right]$
$\displaystyle\quad+\mathcal{S}_{\mathrm{m}}\left[A^{2}(\Psi)\tilde{g}_{\mu\nu};\psi_{\mathrm{m}}\right]\,,$
(14)

where $\Psi$ is the scalar field, $\tilde{R}$ the Ricci scalar,
$A(\Psi)\equiv\varphi^{-1/2}$ the coupling between the scalar field and
matter, and $\mathcal{U}(\Psi)$ the potential Clifton+12

$\mathcal{U}(\Psi)=2\Xi e^{\mp\Psi}\,.$ (15)

The potential $\mathcal{U}(\Psi)=2\Xi e^{-\Psi}$ is the runaway potential of
the chameleon cosmology Clifton+12 ; Joyce+15 ; Khoury2004 . In models endowed
with a chameleon screening, the extra degrees of freedom propagate freely and
mediate a fifth force in regions of weak Newtonian gravitational fields;
however, these degrees of freedom are screened in a high-density environment.
This behaviour coincides with the expected behaviour of the phenomenological
gravitational permittivity of RG. This result therefore suggests a relation
between $\varphi$ and $\epsilon$.

###  III.3 The weak-field limit of CRG

As in standard GR, we take the static, weak-field limit of the metric
$g_{\mu\nu}$,555We adopt the signature $(-,+,+,+)$ for the Minkowski metric.
which can be expanded around the Minkowski metric $\eta_{\mu\nu}$ as

$g_{\mu\nu}\simeq\eta_{\mu\nu}+h_{\mu\nu}\,,$ (16)

where $|h|\ll 1$ and $|\partial_{\mu}h|\ll 1$. We write the metric as

$\displaystyle g_{00}$ $\displaystyle\simeq\eta_{00}-2\Phi\,,$ (17)
$\displaystyle g_{0i}$ $\displaystyle\simeq 0\,,$ (18) $\displaystyle g_{ij}$
$\displaystyle\simeq\eta_{ij}-2U\,,$ (19)

where $\Phi$ and $U$ are two potentials, and we ignore terms of order
$\mathcal{O}(\Phi^{2})\sim\mathcal{O}(U^{2})$. These expressions enable the
calculation of the left-hand side of Eqs. (10) and (11) up to terms of order
$\mathcal{O}(\Phi)\sim\mathcal{O}(U)$ (see Appendix A for the detailed
calculations).

For the right-hand side of Eqs. (10) and (11), we consider a static non-
relativistic fluid: the only non-zero component of the energy-momentum tensor
is $T_{00}\simeq\rho$, and thus its trace is $T\simeq-\rho$. With these
assumptions, the 00-component of the field equations reduces to

$\vec{\nabla}\cdot\left(\varphi\vec{\nabla}\Phi\right)\simeq 8\pi G\rho\,.$
(20)

We thus recover the RG equation if we identify the scalar field with twice the
permittivity: $\varphi=2\epsilon$. In the Newtonian regime, we have a constant
scalar field, namely $\vec{\nabla}\varphi=0$, and thus we recover the standard
Poisson equation for $\varphi=2$.

####  III.3.1 The gravitational field of a spherical source immersed in a
homogeneous background

In Appendix A.1, we compute the gravitational field generated by a spherical
source immersed in a homogeneous background with density $\rho_{\rm bg}$. The
source is described by a density profile $\rho_{s}(r)$ decreasing with $r$. We
estimate the field in the limit $\rho_{s}(r)\gg\rho_{\mathrm{bg}}$, close to
the source, and $\rho_{s}(r)\ll\rho_{\mathrm{bg}}$, at large distance from the
source.

At small distances, we find

$\frac{d\Phi}{dr}=\frac{2}{\varphi}\frac{Gm(<r)}{r^{2}}\,,$ (21)

with the scalar field

$\varphi(r)=2\left[1-\int_{0}^{r}\frac{Gm(<r^{\prime})}{r^{\prime
2}}\,dr^{\prime}\right]\,.$ (22)

At large distances from the source, the scalar field $\varphi$ and the
gravitational field $d\Phi/dr$ satisfy the implicit relation

$\frac{d\ln\varphi}{dr}=\frac{d\Phi}{dr}\left\\{-1-\left[1+\left(\frac{d\Phi}{dr}\right)^{-2}\left(2\Xi-\frac{8\pi
G\rho}{\varphi}\right)\right]^{1/2}\right\\}\,,$ (23)

with $\rho(r)=\rho_{s}(r)+\rho_{\mathrm{bg}}$. This relation sets the
acceleration scale

$a_{\Xi}=\left(2\Xi-\frac{8\pi G\rho}{\varphi}\right)^{1/2}\,.$ (24)

In regions where $d\Phi/dr\gg a_{\Xi}$, the gravitational field
$d\Phi/dr\simeq-(1/2)d\ln\varphi/dr$ has a similar dependence to that at small
distances (Eqs. 21 and 22). In regions where $d\Phi/dr\ll a_{\Xi}$, the RG
acceleration deviates from the Newtonian acceleration. This result resembles
the starting hypothesis of MOND, that introduces the acceleration scale
$a_{0}$ to separate the Newtonian from the modified gravity regimes. Moreover,
in Sect. IV, we find $\Xi\sim\Lambda$. We thus find that $a_{\Xi}\sim
10^{-10}$ m s-2, in the limit $2\Xi\gg 8\pi G\rho/\varphi$, that occurs at
large distances from the source. The existence of this acceleration scale has
been confirmed by a wealth of observations on the scale of galaxies (e.g.
deMartino:2020gfi ; McGaugh:2020ppt ; Chae:2020vlk ; McGaugh:2018 ;
McGaugh:2016leg ). Indeed, the dynamics of disk and elliptical galaxies in the
low-acceleration regime is described in the RG framework without requiring the
existence of dark matter (Cesare+20, ; Cesare+21, ).

The connection between $a_{\Xi}$ and $\Xi$ is similar to the connection
between $a_{0}$ and $\Lambda$ in MOND: $a_{0}\sim\Lambda^{1/2}$. A number of
different sensible arguments have been suggested for the interpretation of the
latter relation Milgrom:1989alt ; Milgrom:1998sy ; Famaey+12 ; Milgrom:2020cch
. In the CRG context, $a_{\Xi}\sim(2\Xi)^{1/2}$ emerges naturally. However,
the physical constant $a_{0}$ introduced by MOND is independent of the source
of the gravitational field, whereas, in CRG, the acceleration scale $a_{\Xi}$
depends on the density field $\rho=\rho_{\mathrm{bg}}+\rho_{s}$ appearing in
Eq. (24): in this limit of large distances from the spherical source,
$\rho_{s}$ acts as a perturbation of $\rho$ and we thus expect a weak
dependence of $a_{\Xi}$ on the source. It remains to be seen whether, when
dropping all the simplifying hypotheses of this section, $a_{\Xi}$ is
consistent with the ubiquitous emergence of $a_{0}$ on the scale of galaxies
in the real Universe.

####  III.3.2 The scalar field $\varphi$ for a spherical source immersed in a
constant background

According to the results of Appendix A.1, the scalar field $\varphi$ is
positive and broadly in the range $[0,\,2]$ consistently with the RG ansatz
$\epsilon\in[0,1]$ for the permittivity. In the limit
$\rho_{\mathrm{bg}}\ll\rho_{s}$, we have

$\varphi_{\mathrm{near}}\simeq
2\left[1-\int_{0}^{r}\frac{Gm(<r^{\prime})}{r^{\prime
2}}\,dr^{\prime}\right]\,$ (25)

whereas in the limit $\rho_{\mathrm{bg}}\gg\rho_{s}$,

$\varphi_{\mathrm{far}}\simeq\frac{4\pi
G}{\Xi}\left(\rho_{\mathrm{bg}}+\rho_{s}\right)\,.$ (26)

These limits are broadly in agreement with the smooth step function considered
for $\epsilon=\epsilon(\rho)$ in refs. Cesare+21 ; Cesare+20 ;
MatsakosDiaferio16 . In those studies, the local mass density was found to be
a good proxy of the transition between the Newtonian and the RG regimes.

###  III.4 Schwarzschild-type solution

In Appendix B, we derive the field equations, Eqs. (116)-(119), for the case
of a static, spherically-symmetric metric of the form

$ds^{2}=-A^{2}dt^{2}+B^{2}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right),$ (27)

where $A=A(r)$ and $B=B(r)$.

Here, we only consider the special case of a point-mass source, so that
$\rho=M\delta^{3}(r)$, within a homogeneous background with density
$\rho_{\mathrm{bg}}$, small enough to provide a negligible contribution to the
gravitational field. We also assume a constant scalar field
$\varphi=\varphi_{0}$, which is a solution of the field equations, Eqs.
(116)-(119). With these assumptions, the field equations simplify to

$\displaystyle\frac{1}{AB}\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}A\right)$
$\displaystyle=\frac{8\pi GM}{\varphi_{0}}\delta^{3}(r)\,,$ (28)
$\displaystyle\vec{\nabla}\cdot\left[\frac{1}{r^{2}}\ln{(AB)}\right]$
$\displaystyle=0\,,$ (29)
$\displaystyle-r\vec{\nabla}\cdot\left(\ln{\frac{A}{B}}\vec{\nabla}\frac{1}{r}\right)+\frac{1}{r^{2}}-\frac{B^{2}}{r^{2}}$
$\displaystyle=0\,,$ (30) $\displaystyle\varphi_{0}$ $\displaystyle=\frac{4\pi
G\rho_{\mathrm{bg}}}{\Xi}\,.$ (31)

Equation (29) yields $B=1/A$, and Eqs. (28) and (29) can then be solved for
$A$:

$A(r)=\left(1-\frac{4GM}{\varphi_{0}r}\right)^{1/2}\,.$ (32)

Similarly to the weak-field limit of a point source described in Sect. III.3
above, we find that the gravitational potential of the source is enhanced by
the factor $2/\varphi_{0}$ (since $0<\varphi_{0}<2$) and the standard GR
solution is recovered when $\varphi_{0}=2$. In the general case of
$\varphi=\varphi(r)$, $A(r)$ will clearly be different, but the role of the
scalar field to enhance the gravitational field, contained in the
$00$-component of the metric, will be preserved. This feature is also evident
from the general Equation (116) in Appendix B: in vacuum this equation reads
$\partial_{r}A=\mathcal{C}B/r^{2}\varphi$, with $\mathcal{C}$ an integration
constant.

##  IV The homogeneous and isotropic Universe in CRG

The covariant formulation of RG enables the description of a homogeneous and
isotropic universe that can be described by the Friedmann-Lemaître-Robertson-
Walker (FLRW) metric

$ds^{2}=-dt^{2}+a(t)^{2}\left(\frac{dr^{2}}{1-kr^{2}}+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta\
d\phi^{2}\right),$ (33)

where $a(t)$ is the scale factor, $k$ the spatial curvature, and
$(t,r,\theta,\phi)$ the co-moving coordinates. By treating the content of the
Universe as a perfect fluid, Eq. (10) readily yields the modified Friedmann
equations (see Appendix C)

$\frac{\ddot{a}}{a}-\frac{1}{3}\left(\frac{\dot{\varphi}^{2}}{\varphi^{2}}-\frac{\ddot{\varphi}}{\varphi}\right)=-\frac{8\pi
G}{3\varphi}\rho$ (34)

and

$\frac{\ddot{a}}{a}+2\frac{\dot{a}^{2}}{a^{2}}+\frac{2k}{a^{2}}+\frac{\dot{a}\dot{\varphi}}{a\varphi}=\frac{8\pi
G}{\varphi}p\,.$ (35)

The equation for the scalar field, Eq. (11), reduces to (see Appendix C)

$\frac{1}{3}\frac{\ddot{\varphi}}{\varphi}+\frac{\dot{a}\dot{\varphi}}{a\varphi}+\frac{2\Xi}{3}=\frac{8\pi
G}{3\varphi}\left(\rho-3p\right)\,.$ (36)

In the JF, the stress-energy tensor is covariantly conserved Faraoni+04 . We
assume the standard equation of state

$p=w\rho\,,$ (37)

with $w=0$ and $w=1/3$ describing the dust and radiation components,
respectively. The dependence of the matter or radiation density on time $t$,
or equivalently on the scale factor $a(t)$, is thus

$\rho(t)=\rho_{0}a^{-3\left(w+1\right)}\,,$ (38)

with $\rho_{0}$ the mean density of the component at present time, $t_{0}$,
when $a(t_{0})=a_{0}=1$.

###  IV.1 A spatially-flat universe: analytic solution

Here, we solve the field equations, Eqs. (34)-(36), in the special case of a
spatially-flat universe with $k=0$. We assume that the universe only contains
baryonic matter with density $\rho_{\rm b}$ and negligible pressure, $p=0$,
namely $w=0$ in the equation of state, Eq. (37).

After substituting Eq. (34) into Eq. (35) and dividing by
$H^{2}\equiv(\dot{a}/a)^{2}$, Eqs. (34)-(36) become

$\displaystyle\frac{1}{H^{2}}\frac{\ddot{a}}{a}-\frac{1}{3H^{2}}\left(\frac{\dot{\varphi}^{2}}{\varphi^{2}}-\frac{\ddot{\varphi}}{\varphi}\right)$
$\displaystyle=-\frac{8\pi G}{3H^{2}\varphi}\rho_{\rm b}\,,$ (39)
$\displaystyle
2+\frac{1}{H}\frac{\dot{\varphi}}{\varphi}+\frac{1}{3H^{2}}\left(\frac{\dot{\varphi}^{2}}{\varphi^{2}}-\frac{\ddot{\varphi}}{\varphi}\right)$
$\displaystyle=\frac{8\pi G}{3H^{2}\varphi}\rho_{\rm b}\,,$ (40)
$\displaystyle\frac{1}{3H^{2}}\frac{\ddot{\varphi}}{\varphi}+\frac{1}{H}\frac{\dot{\varphi}}{\varphi}+\frac{2}{3H^{2}}\Xi$
$\displaystyle=\frac{8\pi G}{3H^{2}\varphi}\rho_{\rm b}\,.$ (41)

The solution of the above system of coupled equations together with the
equation of the mass-energy conservation, Eq. (38), determines the time
evolution of both the scalar field and the scale factor.

We simplify the field equations by introducing the modified cosmological
parameter

$\Omega\equiv\frac{\Omega_{\mathrm{b}}}{\varphi/2}\equiv\frac{16\pi
G}{3H^{2}\varphi}\rho_{\rm b}\,,$ (42)

where $\Omega_{\mathrm{b}}=\rho_{\rm b}/\rho_{\rm cr}\equiv\rho_{\rm b}(8\pi
G/3H^{2})$ is the density parameter associated to the homogeneous baryonic
density $\rho_{\rm b}$. $\Omega$ is analogous to the total density of matter
in the standard model, that includes both baryonic and dark matter; in CRG,
the gravitational role of the baryonic matter density is amplified by the
factor $2/\varphi$, unlike the standard model, where the baryonic and the dark
matter components simply add up.

Similarly, we define

$\Omega_{\Xi}=\frac{\Xi}{3H^{2}}\,.$ (43)

$\Omega_{\Xi}$ is analogous to the standard cosmological parameter
$\Omega_{\Lambda}=\Lambda/3H^{2}$, with $\Lambda$ the cosmological constant in
the GR field equations; therefore, in CRG, $\Xi$ exactly plays the role of
$\Lambda$ in the standard model. At $t=t_{0}$, the values of the two
parameters $\Omega$ and $\Omega_{\Xi}$ are $\Omega_{0}=16\pi G\rho_{\rm
b0}/3H_{0}^{2}\varphi_{0}$ and $\Omega_{\Xi 0}=\Xi/3H_{0}^{2}$.

We further introduce two deceleration parameters, $q$ and $q_{\varphi}$,
related to the scale factor and the scalar field, respectively:

$q\equiv-\frac{a\ddot{a}}{\dot{a}^{2}}=-\frac{\ddot{a}}{aH^{2}}\,,\quad
q_{\varphi}\equiv-\frac{\ddot{\varphi}}{\varphi H^{2}}\,.$ (44)

With the above definitions, the field equations, Eqs. (39)-(41), become

$\displaystyle-q-\frac{1}{3H^{2}}\left(\frac{\dot{\varphi}}{\varphi}\right)^{2}-\frac{1}{3}q_{\varphi}$
$\displaystyle=-\frac{1}{2}\Omega\,,$ (45) $\displaystyle
2+\frac{1}{H}\frac{\dot{\varphi}}{\varphi}+\frac{1}{3H^{2}}\left(\frac{\dot{\varphi}}{\varphi}\right)^{2}+\frac{1}{3}q_{\varphi}$
$\displaystyle=\frac{1}{2}\Omega\,,$ (46)
$\displaystyle-\frac{1}{3}q_{\varphi}+\frac{1}{H}\frac{\dot{\varphi}}{\varphi}+2\Omega_{\Xi}$
$\displaystyle=\frac{1}{2}\Omega\,.$ (47)

It is convenient to rewrite $\dot{\varphi}/H\varphi$ in terms of the new
quantity

$\zeta\equiv\frac{\dot{\varphi}}{H\varphi}\,.$ (48)

The field equations, Eqs. (45)-(47), now become

$\displaystyle q+\frac{1}{3}\zeta^{2}+\frac{1}{3}q_{\varphi}$
$\displaystyle=\frac{1}{2}\Omega\,,$ (49) $\displaystyle
2+\zeta+\frac{1}{3}\zeta^{2}+\frac{1}{3}q_{\varphi}$
$\displaystyle=\frac{1}{2}\Omega\,,$ (50)
$\displaystyle-\frac{1}{3}q_{\varphi}+\zeta$
$\displaystyle=\frac{1}{2}\Omega-2\Omega_{\Xi}\,.$ (51)

Combining the first two equations yields

$q-\zeta=2\,,$ (52)

namely, with the definitions of $q$ and $\zeta$,

$\frac{d}{dt}\ln\varphi=-\frac{d}{dt}\left(\ln\dot{a}+2\ln{a}\right)\,.$ (53)

This equation can be integrated to obtain
$\varphi=\mathcal{C}/(a^{2}\dot{a})$, with $\mathcal{C}$ an integration
constant. Its value can be found from the boundary conditions $a(t_{0})=1$,
$H(t_{0})=H_{0}$, and $\varphi(t_{0})=\varphi_{0}$. We obtain
$\mathcal{C}=H_{0}\varphi_{0}$ and thus

$\varphi=\frac{H_{0}\varphi_{0}}{Ha^{3}}\,.$ (54)

Moreover, since $3Ha^{3}=da^{3}/dt$ we can write the scalar field as

$\varphi=3H_{0}\varphi_{0}\left(\frac{da^{3}}{dt}\right)^{-1}\,.$ (55)

This expression shows that the scalar field $\varphi$ is inversely
proportional to the rate of the variation of the volume of the universe.

As shown in Appendix D, adding up Eqs. (49) and (51) yields two solutions for
$dH/dt$

$\frac{dH}{dt}=\pm\sqrt{3}\left(H^{4}+\Omega_{0}H_{0}H^{3}-2\Omega_{\Xi
0}H_{0}^{2}H^{2}\right)^{1/2}\,,$ (56)

each being a first order differential equation for $H$. Their integration
yields

$H(t)=\frac{4\Omega_{\Xi
0}H_{0}}{\Omega_{0}\mp\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}\sin\left[\sqrt{6\Omega_{\Xi 0}}(H_{0}t+C_{1})\right]}\,,$ (57)

with the upper and lower signs corresponding to those of Eq. (56).666This
convention holds throughout the remaining of this section. The integration
constant, $C_{1}$, can be found from the boundary conditions at $t=t_{0}$:

$C_{1}=-H_{0}t_{0}\pm\frac{1}{\sqrt{6\Omega_{\Xi
0}}}\arcsin\left[\frac{\Omega_{0}-4\Omega_{\Xi
0}}{\left(\Omega_{0}^{2}+8\Omega_{\Xi 0}\right)^{1/2}}\right]\,.$ (58)

Furthermore, the scale factor can be derived by integrating Eq. (57)

$\displaystyle a(t)=$ $\displaystyle
C_{2}\left(\left|\frac{\Omega_{0}\tan\tau\mp\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}-(8\Omega_{\Xi
0})^{1/2}}{\Omega_{0}\tan\tau\mp\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}+(8\Omega_{\Xi 0})^{1/2}}\right|\right)^{1/\sqrt{3}}\,,$ (59)

where

$\tau=\left(\frac{3}{2}\Omega_{\Xi 0}\right)^{1/2}(H_{0}t+C_{1})\,,$ (60)

and $C_{2}$ an integration constant, whose value can be fixed when $a_{0}=1$:

$C_{2}=\left(\left|\frac{\Omega_{0}\tan\tau_{0}\mp\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}+(8\Omega_{\Xi
0})^{1/2}}{\Omega_{0}\tan\tau_{0}\mp\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}-(8\Omega_{\Xi 0})^{1/2}}\right|\right)^{1/\sqrt{3}}\,,$ (61)

with

$\tau_{0}=\pm\frac{1}{2}\arcsin\left[\frac{\Omega_{0}-4\Omega_{\Xi
0}}{\left(\Omega_{0}^{2}+8\Omega_{\Xi 0}\right)^{1/2}}\right]\,.$ (62)

Inverting Eq. (62) yields

$\sin\left(\pm 2\tau_{0}\right)=\frac{\Omega_{0}-4\Omega_{\Xi
0}}{\left(\Omega_{0}^{2}+8\Omega_{\Xi 0}\right)^{1/2}},$ (63)

which implies that the cosmological parameters $\Omega$ and $\Omega_{\Xi}$
satisfy the relation

$-1\leq\frac{\Omega_{0}-4\Omega_{\Xi 0}}{\left(\Omega_{0}^{2}+8\Omega_{\Xi
0}\right)^{1/2}}\leq 1.$ (64)

$\Omega_{0}$ must be positive, because it is proportional to the energy
density $\rho_{\rm b0}$ according to Eq. (42); therefore, the only physically
viable solutions of these inequalities are

$\Omega_{0}>0;\qquad 0<\Omega_{\Xi 0}\leq\frac{\Omega_{0}+1}{2}\,.$ (65)

The parameter $\Xi$ must thus be positive.

####  IV.1.1 The supernovae data in the CRG universe

The first test of the viability of CRG is its ability to properly describe the
Hubble diagram of the observed Type Ia Supernovae (SNIa) at high redshift.
Using the scale factor in Eq. (59), the luminosity distance can be computed
exactly, for both signs appearing in Eq. (59) (see Appendix E). The upper sign
is referred to as the $CRG+$ solution, while the lower as the $CRG-$ solution.

As in the standard cosmological model, the values of the cosmological
parameters can be set by fitting the SNIa Hubble diagram. Here, we limit our
analysis to simply show that a set of viable values of these parameters is
consistent with the observations. We postpone a rigorous determination of
their values to future work.

The CRG parameters defined in Eqs. (42)-(43) play roughly the same role as the
parameters $\Omega_{0}$ and $\Omega_{\Lambda 0}$ of the standard model. To
test the ability of CRG to reproduce the Supernovae data, we can thus simply
adopt the corresponding values derived in the standard model. A number of
observations, including data from the Planck mission Aghanim:2018eyx , suggest
$\Omega_{0}=0.31$. According to the limits of Eq. (65), we set $\Omega_{\Xi
0}=0.65$. For the value of the Hubble constant we adopt $H_{0}=67.7$ km s-1
Mpc-1 Aghanim:2018eyx .

Figure 1 compares the distance modulus $m-M$ expected in CRG, according to the
luminosity distance of Eq. (154) derived with this combination of parameters,
with the SNIa data, taken from the Supernova Cosmology Project Union 2.1
Compilation Amanullah2010 . The qualitative agreement between the data and the
CRG distance modulus suggests that CRG can indeed describe the data for a
specific combination of the values of $\Omega_{0}$, $\Omega_{\Xi 0}$ and
$H_{0}$. The agreement shown in Fig. 1 also justifies the adoption of the
value of $\Xi\sim\Lambda$ when, in Sect. III.3.1, we estimated the
acceleration scale $(2\Xi)^{1/2}\sim 10^{-10}$ m s-2 that separates the
Newtonian and the RG regime in the weak-field limit.

Figure 1: Distance modulus computed with the ”CRG+” solution (solid green
line) and the ”CRG-” solution (dashed red line) compared with supernovae data
from the Supernova Cosmology Project Union 2.1 Compilation Amanullah2010 . For
the two curves, we adopt $H_{0}=67.7$ km s-1 Mpc-1, $\Omega_{\Xi 0}=0.65$, and
$\Omega_{0}=0.31$.

##  V Conclusion

RG was originally introduced to describe the dynamics of galaxies and galaxy
systems without the aid of dark matter. The contribution of dark matter to the
gravitational field is mimicked by the gravitational permittivity
$\epsilon(\rho)$, a monotonic function of the mass density $\rho$
MatsakosDiaferio16 . Here, we propose a covariant extension of RG, CRG, and
show that CRG belongs to the family of scalar-tensor theories, thus inheriting
all their general properties. A scalar-tensor theory is specified by the self-
interacting potential $\mathcal{V}(\varphi)$ of the scalar field $\varphi$ and
by the general differentiable function $\mathcal{W}(\varphi)$ appearing in the
Lagrangian density. For CRG, we adopt $\mathcal{V}(\varphi)=-\Xi\varphi$, with
$\Xi$ a normalization constant, and $\mathcal{W}(\varphi)=-1$. In the weak-
field limit, this theory correctly reduces to the original phenomenological RG
and identifies the gravitational permittivity with the scalar field, with
$\varphi=2\epsilon$.

When used to describe an expanding and homogeneous universe, the scalar field
$\varphi$ is also responsible for the observed accelerated expansion of the
Universe. The cosmological density associated to the cosmological constant in
the standard model, $\Omega_{\Lambda 0}=\Lambda/3H_{0}^{2}$, is now replaced
by the CRG parameter $\Omega_{\Xi 0}=\Xi/3H_{0}^{2}$. It follows that the
normalization constant $\Xi$ plays the role of the cosmological constant of
the standard model and might thus also assume a value comparable to that of
$\Lambda$. In addition, in the weak-field limit, $\Xi$ sets an acceleration
scale, $(2\Xi)^{1/2}\sim 10^{-10}$ m s-2, below which RG deviates from
Newtonian gravity and appears to describe the dynamics of disk and elliptical
galaxies without the aid of dark matter Cesare+20 ; Cesare+21 . This
acceleration scale is indeed present in real systems Chae:2020vlk ;
McGaugh:2018  and is comparable to the acceleration $a_{0}$ introduced in MOND
Milgrom1 ; Milgrom2 ; Milgrom3 . Therefore, CRG provides a connection between
phenomena generally attributed to dark matter and dark energy separately. In
addition, the known relation $a_{0}\sim\Lambda^{1/2}$ Milgrom:1998sy ;
Milgrom:2020cch  appears naturally.

The feature of connecting dark matter and dark energy within a unified
scenario is not unique to CRG. The same properties are currently shared by
other models, including some $f(R)$ theories f(R)Faraoni , quartessence
theories Brandenberger+19  or generalized Chaplygin gas models Bento:2002ps ;
Zhang:2004gc .

Here, we showed that CRG can correctly reproduce the Hubble diagram of SNIa.
Additional probes are clearly necessary to assess if the evolution of the
universe described by the present model is consistent with observations. The
power spectrum of the temperature anisotropies of the Cosmic Microwave
Background and the power spectrum of matter density perturbations are required
to constrain the Lagrangian density and its parameters (see, e.g.,
Noller:2018wyv ; Huterer:2013xky ). Similarly, the analysis of the evolution
of the density perturbations, at least in the linear perturbation theory
DiPorto:2007ovd ; BuenoSanchez:2011zz ; Kofinas:2017clq , and of the role of
the scalar field and its perturbation on structure formation are of crucial
importance to test the viability of CRG.

In addition, since the scalar field drives the accelerated expansion of the
Universe, we need to investigate the connection of CRG with current dark
energy models Pettorino:2005pv ; Capozziello:2005mj ; Frusciante:2019xia .
This study might provide specific predictions of the redshift evolution of an
effective dark energy fluid in CRG. When tested with measures from the
upcoming Euclid mission Amendola:2016saw , these predictions could
discriminate CRG from current dark energy models. We plan to tackle all these
issues in future work.

###### Acknowledgements.

We thank Mariano Cadoni for useful discussions. An early version of this work
was the Master thesis in Physics at the University of Torino of APS. We
acknowledge partial support from the INFN grant InDark and the Italian
Ministry of Education, University and Research (MIUR) under the Departments of
Excellence grant L.232/2016. This research has made use of NASA’s Astrophysics
Data System Bibliographic Services.

## References

  * (1) J. P. Ostriker and P. J. Steinhardt, “The observational case for a low-density Universe with a non-zero cosmological constant,” Nature (London), vol. 377, pp. 600–602, Oct. 1995.
  * (2) R. H. Cyburt, B. D. Fields, K. A. Olive, and T.-H. Yeh, “Big bang nucleosynthesis: Present status,” Reviews of Modern Physics, vol. 88, p. 015004, Jan. 2016.
  * (3) N. Aghanim et al., “Planck 2018 results. V. CMB power spectra and likelihoods,” Astron. Astrophys., vol. 641, p. A5, 2020.
  * (4) D. Clowe, M. Bradac, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, “A direct empirical proof of the existence of dark matter,” Astrophys. J. Lett., vol. 648, pp. L109–L113, 2006.
  * (5) S. Dodelson et al., “The Three-Dimensional Power Spectrum from Angular Clustering of Galaxies in early SDSS Data,” Astrophys. J., vol. 572, pp. 140–156, 2001.
  * (6) A. Del Popolo, “Nonbaryonic Dark Matter in Cosmology,” Int. J. Mod. Phys. D, vol. 23, p. 1430005, 2014.
  * (7) Y. Akrami et al., “Planck 2018 results. IV. Diffuse component separation,” Astron. Astrophys., vol. 641, p. A4, 2020.
  * (8) N. Aghanim et al., “Planck 2018 results. I. Overview and the cosmological legacy of Planck,” Astron. Astrophys., vol. 641, p. A1, 2020\.
  * (9) A. G. Riess et al., “Observational evidence from supernovae for an accelerating universe and a cosmological constant,” Astron. J., vol. 116, pp. 1009–1038, 1998.
  * (10) S. Perlmutter et al., “Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae,” Astrophys. J., vol. 517, pp. 565–586, 1999\.
  * (11) N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys., vol. 641, p. A6, 2020.
  * (12) A. G. Riess, L. M. Macri, S. L. Hoffmann, D. Scolnic, S. Casertano, A. V. Filippenko, B. E. Tucker, M. J. Reid, D. O. Jones, J. M. Silverman, R. Chornock, P. Challis, W. Yuan, P. J. Brown, and R. J. Foley, “A 2.4% Determination of the Local Value of the Hubble Constant,” Astrophys. J. , vol. 826, p. 56, July 2016.
  * (13) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, “Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM,” Astrophys. J. , vol. 876, p. 85, May 2019.
  * (14) L. Verde, T. Treu, and A. Riess, “Tensions between the Early and the Late Universe,” 7 2019.
  * (15) G. Mathews, A. Kedia, N. Sasankan, M. Kusakabe, Y. Luo, T. Kajino, D. Yamazaki, T. Makki, and M. Eid, “Cosmological Solutions to the Lithium Problem,” JPS Conf. Proc., vol. 31, p. 011033, 2020.
  * (16) P. Ade et al., “Planck 2015 results. XVI. Isotropy and statistics of the CMB,” Astron. Astrophys., vol. 594, p. A16, 2016.
  * (17) D. J. Schwarz, C. J. Copi, D. Huterer, and G. D. Starkman, “CMB Anomalies after Planck,” Class. Quant. Grav., vol. 33, no. 18, p. 184001, 2016.
  * (18) E. Di Valentino, A. Melchiorri, and J. Silk, “Planck evidence for a closed Universe and a possible crisis for cosmology,” Nature Astron., vol. 4, no. 2, pp. 196–203, 2019.
  * (19) M. Raveri, “Are cosmological data sets consistent with each other within the $\Lambda$ cold dark matter model?,” Phys. Rev. D, vol. 93, no. 4, p. 043522, 2016.
  * (20) H. Hildebrandt et al., “KiDS-450: Cosmological parameter constraints from tomographic weak gravitational lensing,” Mon. Not. Roy. Astron. Soc., vol. 465, p. 1454, 2017.
  * (21) S. Weinberg, “The Cosmological Constant Problem,” Rev. Mod. Phys., vol. 61, pp. 1–23, 1989.
  * (22) A. Padilla, “Lectures on the Cosmological Constant Problem,” 2 2015.
  * (23) A. Joyce, B. Jain, J. Khoury, and M. Trodden, “Beyond the Cosmological Standard Model,” Phys. Rept., vol. 568, pp. 1–98, 2015.
  * (24) P. Peebles and B. Ratra, “The Cosmological Constant and Dark Energy,” Rev. Mod. Phys., vol. 75, pp. 559–606, 2003.
  * (25) E. J. Copeland, M. Sami, and S. Tsujikawa, “Dynamics of dark energy,” Int. J. Mod. Phys. D, vol. 15, pp. 1753–1936, 2006.
  * (26) L. Amendola et al., “Cosmology and fundamental physics with the Euclid satellite,” Living Rev. Rel., vol. 21, no. 1, p. 2, 2018.
  * (27) A. Del Popolo and M. Le Delliou, “Small scale problems of the $\Lambda$CDM model: a short review,” Galaxies, vol. 5, no. 1, p. 17, 2017.
  * (28) P. Kroupa, “The dark matter crisis: falsification of the current standard model of cosmology,” Publ. Astron. Soc. Austral., vol. 29, pp. 395–433, 2012.
  * (29) I. de Martino, S. S. Chakrabarty, V. Cesare, A. Gallo, L. Ostorero, and A. Diaferio, “Dark Matters on the Scale of Galaxies,” Universe, vol. 6, p. 107, Aug. 2020.
  * (30) G. Bertone, D. Hooper, and J. Silk, “Particle dark matter: Evidence, candidates and constraints,” Phys. Rept., vol. 405, pp. 279–390, 2005\.
  * (31) M. Tanabashi et al., “Review of Particle Physics,” Phys. Rev. D, vol. 98, no. 3, p. 030001, 2018.
  * (32) M. R. Buckley and A. H. G. Peter, “Gravitational probes of dark matter physics,” Phys. Rept., vol. 761, pp. 1–60, 2018.
  * (33) L. Berezhiani and J. Khoury, “Theory of dark matter superfluidity,” Phys. Rev. D, vol. 92, p. 103510, 2015.
  * (34) L. Blanchet, “Gravitational polarization and the phenomenology of MOND,” Class. Quant. Grav., vol. 24, pp. 3529–3540, 2007.
  * (35) L. Blanchet and L. Heisenberg, “Dipolar Dark Matter as an Effective Field Theory,” Phys. Rev. D, vol. 96, no. 8, p. 083512, 2017.
  * (36) M. Milgrom, “A Modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis,” Astrophys. J., vol. 270, pp. 365–370, 1983.
  * (37) M. Milgrom, “A Modification of the Newtonian dynamics: Implications for galaxies,” Astrophys. J., vol. 270, pp. 371–383, 1983.
  * (38) M. Milgrom, “A modification of the Newtonian dynamics: implications for galaxy systems,” Astrophys. J., vol. 270, pp. 384–389, 1983.
  * (39) R. H. Sanders and S. S. McGaugh, “Modified Newtonian dynamics as an alternative to dark matter,” Ann. Rev. Astron. Astrophys., vol. 40, pp. 263–317, 2002.
  * (40) B. Famaey and S. McGaugh, “Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions,” Living Rev. Rel., vol. 15, p. 10, 2012.
  * (41) M. Milgrom, “MOND theory,” Can. J. Phys., vol. 93, no. 2, pp. 107–118, 2015.
  * (42) C. Skordis and T. Złośnik, “Gravitational alternatives to dark matter with tensor mode speed equaling the speed of light,” Phys. Rev. D, vol. 100, no. 10, p. 104013, 2019.
  * (43) T. G. Złośnik and C. Skordis, “Cosmology of the Galileon extension of Bekenstein’s theory of relativistic modified Newtonian dynamics,” Phys. Rev. D, vol. 95, no. 12, p. 124023, 2017.
  * (44) E. G. Ferreira, G. Franzmann, J. Khoury, and R. Brandenberger, “Unified Superfluid Dark Sector,” JCAP, vol. 08, p. 027, 2019.
  * (45) A. Sakharov, “Vacuum quantum fluctuations in curved space and the theory of gravitation,” Usp. Fiz. Nauk, vol. 161, no. 5, pp. 64–66, 1991.
  * (46) T. Padmanabhan, “Emergent Gravity Paradigm: Recent Progress,” Mod. Phys. Lett. A, vol. 30, no. 03n04, p. 1540007, 2015.
  * (47) E. P. Verlinde, “Emergent Gravity and the Dark Universe,” SciPost Phys., vol. 2, no. 3, p. 016, 2017.
  * (48) M. Cadoni, R. Casadio, A. Giusti, W. Mück, and M. Tuveri, “Effective Fluid Description of the Dark Universe,” Phys. Lett. B, vol. 776, pp. 242–248, 2018.
  * (49) M. Cadoni, R. Casadio, A. Giusti, and M. Tuveri, “Emergence of a Dark Force in Corpuscular Gravity,” Phys. Rev. D, vol. 97, no. 4, p. 044047, 2018\.
  * (50) M. Tuveri and M. Cadoni, “Galactic dynamics and long-range quantum gravity,” Phys. Rev. D, vol. 100, no. 2, p. 024029, 2019.
  * (51) M. Cadoni, A. P. Sanna, and M. Tuveri, “Anisotropic fluid cosmology: An alternative to dark matter?,” Phys. Rev. D, vol. 102, no. 2, p. 023514, 2020.
  * (52) T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, “Modified gravity and cosmology,” Physics Reports, vol. 513, pp. 1–189, Mar 2012\.
  * (53) D. Lovelock, “The Einstein tensor and its generalizations,” J. Math. Phys., vol. 12, pp. 498–501, 1971.
  * (54) D. Lovelock, “The four-dimensionality of space and the einstein tensor,” J. Math. Phys., vol. 13, pp. 874–876, 1972.
  * (55) P. D. Mannheim and D. Kazanas, “Newtonian limit of conformal gravity and the lack of necessity of the second order Poisson equation,” Gen. Rel. Grav., vol. 26, pp. 337–361, 1994.
  * (56) C. M. Bender and P. D. Mannheim, “No-ghost theorem for the fourth-order derivative Pais-Uhlenbeck oscillator model,” Phys. Rev. Lett., vol. 100, p. 110402, 2008.
  * (57) P. D. Mannheim, “Making the Case for Conformal Gravity,” Found. Phys., vol. 42, pp. 388–420, 2012.
  * (58) A. Diaferio, L. Ostorero, and V. F. Cardone, “Gamma-ray bursts as cosmological probes: LambdaCDM vs. conformal gravity,” JCAP, vol. 10, p. 008, 2011.
  * (59) D. Elizondo and G. Yepes, “Can conformal Weyl gravity be considered a viable cosmological theory?,” Astrophys. J., vol. 428, pp. 17–20, 1994.
  * (60) K. Horne, “X-ray gas in the galaxy cluster abell 2029: conformal gravity versus dark matter,” Monthly Notices of the Royal Astronomical Society, vol. 369, no. 4, pp. 1667–1676, 2006.
  * (61) A. Diaferio and L. Ostorero, “X-ray clusters of galaxies in conformal gravity,” Mon. Not. Roy. Astron. Soc., vol. 393, p. 215, 2009.
  * (62) M. C. Campigotto, A. Diaferio, and L. Fatibene, “Conformal gravity: light deflection revisited and the galactic rotation curve failure,” Classical and Quantum Gravity, vol. 36, p. 245014, Dec. 2019.
  * (63) S. Carneiro, “Can the cosmological dark sector be modeled by a single scalar field?,” Gen. Rel. Grav., vol. 50, no. 9, p. 114, 2018.
  * (64) D. Bertacca, N. Bartolo, and S. Matarrese, “Unified Dark Matter Scalar Field Models,” Adv. Astron., vol. 2010, p. 904379, 2010.
  * (65) S. Camera, T. D. Kitching, A. F. Heavens, D. Bertacca, and A. Diaferio, “Measuring Unified Dark Matter with 3D cosmic shear,” Mon. Not. Roy. Astron. Soc., vol. 415, pp. 399–409, 2011.
  * (66) S. Camera, M. Martinelli, and D. Bertacca, “Does quartessence ease cosmic tensions?,” Phys. Dark Univ., vol. 23, p. 100247, 2019.
  * (67) D. Bertacca, N. Bartolo, A. Diaferio, and S. Matarrese, “How the Scalar Field of Unified Dark Matter Models Can Cluster,” JCAP, vol. 10, p. 023, 2008\.
  * (68) S. Camera, D. Bertacca, A. Diaferio, N. Bartolo, and S. Matarrese, “Weak lensing signal in Unified Dark Matter models,” Mon. Not. Roy. Astron. Soc., vol. 399, pp. 1995–2003, 2009.
  * (69) Y. Fujii and K. Maeda, The scalar-tensor theory of gravitation.  Cambridge Monographs on Mathematical Physics, Cambridge University Press, 7 2007.
  * (70) I. Quiros, “Selected topics in scalar-tensor theories and beyond,” International Journal of Modern Physics D, vol. 28, pp. 1930012–156, Jan. 2019\.
  * (71) T. Matsakos and A. Diaferio, “Dynamics of galaxies and clusters in refracted gravity,” arXiv:1603.04943 [astro-ph.GA].
  * (72) V. Cesare, A. Diaferio, T. Matsakos, and G. Angus, “Dynamics of DiskMass Survey galaxies in refracted gravity,” Astron. Astrophys., vol. 637, p. A70, 2020.
  * (73) M. A. Bershady, M. A. Verheijen, R. A. Swaters, D. R. Andersen, K. B. Westfall, and T. Martinsson, “The DiskMass Survey. I. Overview,” Astrophys. J., vol. 716, pp. 198–233, 2010.
  * (74) V. Cesare, A. Diaferio, and T. Matsakos, “The dynamics of three nearby E0 galaxies in refracted gravity,” arXiv e-prints, p. arXiv:2102.12499, Feb. 2021.
  * (75) V. Faraoni, Cosmology in scalar tensor gravity, vol. 139.  2004\.
  * (76) S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity.  New York: John Wiley and Sons, 1972.
  * (77) I. Quiros, R. García-Salcedo, T. Gonzalez, F. Horta-Rangel, and J. Saavedra, “Brans-Dicke Galileon and the variational principle,” Eur. J. Phys., vol. 37, no. 5, p. 055605, 2016.
  * (78) Y.-C. Li, F.-Q. Wu, and X. Chen, “Constraints on the Brans-Dicke gravity theory with the Planck data,” Phys. Rev. D, vol. 88, p. 084053, 2013.
  * (79) A. Avilez and C. Skordis, “Cosmological constraints on Brans-Dicke theory,” Phys. Rev. Lett., vol. 113, no. 1, p. 011101, 2014.
  * (80) O. Hrycyna, M. Szydłowski, and M. Kamionka, “Dynamics and cosmological constraints on Brans-Dicke cosmology,” Phys. Rev. D, vol. 90, no. 12, p. 124040, 2014.
  * (81) T. Kobayashi, “Horndeski theory and beyond: a review,” Rept. Prog. Phys., vol. 82, no. 8, p. 086901, 2019.
  * (82) R. Kase and S. Tsujikawa, “Dark energy in Horndeski theories after GW170817: A review,” Int. J. Mod. Phys. D, vol. 28, no. 05, p. 1942005, 2019.
  * (83) H. Wang et al., “The GW170817/GRB 170817A/AT 2017gfo Association: Some Implications for Physics and Astrophysics,” Astrophys. J. Lett., vol. 851, no. 1, p. L18, 2017.
  * (84) J. Sakstein and B. Jain, “Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories,” Phys. Rev. Lett., vol. 119, no. 25, p. 251303, 2017.
  * (85) J. Noller and A. Nicola, “Cosmological parameter constraints for Horndeski scalar-tensor gravity,” Phys. Rev. D, vol. 99, no. 10, p. 103502, 2019\.
  * (86) P. Creminelli and F. Vernizzi, “Dark Energy after GW170817 and GRB170817A,” Phys. Rev. Lett., vol. 119, no. 25, p. 251302, 2017.
  * (87) J. Sakstein and B. Jain, “Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories,” Phys. Rev. Lett., vol. 119, no. 25, p. 251303, 2017.
  * (88) T. P. Sotiriou and V. Faraoni, “f(R) Theories Of Gravity,” Rev. Mod. Phys., vol. 82, pp. 451–497, 2010.
  * (89) J. Gleyzes, D. Langlois, F. Piazza, and F. Vernizzi, “Healthy theories beyond Horndeski,” Phys. Rev. Lett., vol. 114, no. 21, p. 211101, 2015.
  * (90) J. Khoury and A. Weltman, “Chameleon cosmology,” Phys. Rev. D, vol. 69, p. 044026, 2004.
  * (91) S. McGaugh, “Predictions and Outcomes for the Dynamics of Rotating Galaxies,” Galaxies, vol. 8, no. 2, p. 35, 2020.
  * (92) K.-H. Chae, M. Bernardi, H. Domínguez Sánchez, and R. K. Sheth, “On the Presence of a Universal Acceleration Scale in Elliptical Galaxies,” Astrophys. J. Lett., vol. 903, no. 2, p. L31, 2020.
  * (93) S. S. McGaugh, P. Li, F. Lelli, and J. M. Schombert, “Presence of a fundamental acceleration scale in galaxies,” Nature Astronomy, vol. 2, pp. 924–924, Nov. 2018.
  * (94) S. McGaugh, F. Lelli, and J. Schombert, “Radial Acceleration Relation in Rotationally Supported Galaxies,” Phys. Rev. Lett., vol. 117, no. 20, p. 201101, 2016.
  * (95) M. Milgrom, “Alternatives to Dark Matter,” 2 1989.
  * (96) M. Milgrom, “The modified dynamics as a vacuum effect,” Phys. Lett. A, vol. 253, pp. 273–279, 1999.
  * (97) M. Milgrom, “The $a_{0}$ – cosmology connection in MOND,” arXiv:2001.09729 [astro-ph.GA].
  * (98) R. Amanullah et al., “Spectra and Light Curves of Six Type Ia Supernovae at 0.511 ¡ z ¡ 1.12 and the Union2 Compilation,” Astrophys. J., vol. 716, pp. 712–738, 2010.
  * (99) R. Brandenberger, R. R. Cuzinatto, J. Fröhlich, and R. Namba, “New Scalar Field Quartessence,” JCAP, vol. 02, p. 043, 2019.
  * (100) M. Bento, O. Bertolami, and A. Sen, “Generalized Chaplygin gas, accelerated expansion and dark energy matter unification,” Phys. Rev. D, vol. 66, p. 043507, 2002.
  * (101) X. Zhang, F.-Q. Wu, and J. Zhang, “A New generalized Chaplygin gas as a scheme for unification of dark energy and dark matter,” JCAP, vol. 01, p. 003, 2006.
  * (102) J. Noller and A. Nicola, “Cosmological parameter constraints for Horndeski scalar-tensor gravity,” Phys. Rev. D, vol. 99, no. 10, p. 103502, 2019\.
  * (103) D. Huterer, D. Kirkby, R. Bean, A. Connolly, K. Dawson, S. Dodelson, A. Evrard, B. Jain, M. Jarvis, E. Linder, R. Mandelbaum, M. May, A. Raccanelli, B. Reid, E. Rozo, F. Schmidt, N. Sehgal, A. Slosar, A. van Engelen, H.-Y. Wu, and G. Zhao, “Growth of cosmic structure: Probing dark energy beyond expansion,” Astroparticle Physics, vol. 63, pp. 23–41, Mar. 2015.
  * (104) C. Di Porto and L. Amendola, “Observational constraints on the linear fluctuation growth rate,” Phys. Rev. D, vol. 77, p. 083508, 2008.
  * (105) J. Bueno Sanchez and L. Perivolaropoulos, “Dark energy and matter perturbations in scalar-tensor theories of gravity,” J. Phys. Conf. Ser., vol. 283, p. 012006, 2011.
  * (106) G. Kofinas and N. A. Lima, “Dynamics of cosmological perturbations in modified Brans-Dicke cosmology with matter-scalar field interaction,” Phys. Rev. D, vol. 96, no. 8, p. 084016, 2017.
  * (107) V. Pettorino, C. Baccigalupi, and F. Perrotta, “Scaling solutions in scalar-tensor cosmologies,” JCAP, vol. 12, p. 003, 2005.
  * (108) S. Capozziello, S. Nojiri, and S. Odintsov, “Dark energy: The Equation of state description versus scalar-tensor or modified gravity,” Phys. Lett. B, vol. 634, pp. 93–100, 2006.
  * (109) N. Frusciante and L. Perenon, “Effective field theory of dark energy: A review,” Phys. Rept., vol. 857, pp. 1–63, 2020.

##  Appendix A Derivation of the weak-field limit

Adopting the static, weak-field metric components of Eqs. (17)-(19), the
Christoffel symbols and the resulting components of the Ricci tensor are

$\displaystyle\Gamma_{00}^{0}$
$\displaystyle=\Gamma_{ij}^{0}=\Gamma_{0j}^{i}=0\,,$ (66)
$\displaystyle\Gamma_{0i}^{0}$ $\displaystyle\simeq\partial_{i}\Phi\,,$ (67)
$\displaystyle\Gamma_{00}^{i}$
$\displaystyle\simeq\delta^{ij}\partial_{j}\Phi\,,$ (68)
$\displaystyle\Gamma_{jk}^{i}$
$\displaystyle\simeq-\left(\delta_{k}^{i}\partial_{j}U+\delta_{j}^{i}\partial_{k}U-\delta^{il}\delta_{jk}\partial_{l}U\right)\,,$
(69) $\displaystyle R_{00}$ $\displaystyle\simeq-\nabla^{2}\Phi\,,$ (70)
$\displaystyle R_{0i}$ $\displaystyle=0\,,$ (71) $\displaystyle R_{ij}$
$\displaystyle\simeq\partial_{i}\partial_{j}(\Phi-U)-\delta_{ij}\nabla^{2}U\,.$
(72)

The $\nabla_{\mu}\nabla_{\nu}\varphi$ components are

$\displaystyle\nabla_{0}\nabla_{0}\varphi\simeq-\vec{\nabla}\Phi\cdot\vec{\nabla}\varphi,$
(73) $\displaystyle\nabla_{0}\nabla_{i}\varphi=0,$ (74)
$\displaystyle\nabla_{i}\nabla_{j}\varphi\simeq\partial_{i}\partial_{j}\varphi+\partial_{j}U\partial_{i}\varphi+\partial_{i}U\partial_{j}\varphi$
$\displaystyle\qquad\qquad-\delta_{ij}\vec{\nabla}U\cdot\vec{\nabla}\varphi\,,$
(75)

and

$\displaystyle\nabla^{\alpha}\nabla_{\alpha}\varphi\simeq\vec{\nabla}(\Phi-U)\cdot\vec{\nabla}\varphi+\nabla^{2}\varphi\,.$
(76)

For a static non-relativistic fluid, whose pressure $p$ is negligible compared
to its density $\rho$, $p\ll\rho$, the equation $u^{\alpha}u_{\alpha}=-1$
implies $u_{0}=(-g_{00})^{1/2}$; therefore, the components of the energy-
momentum tensor $T_{\mu\nu}=\rho u_{\mu}u_{\nu}$ are

$\displaystyle T_{00}$ $\displaystyle\simeq\rho,$ (77) $\displaystyle T_{0i}$
$\displaystyle=0,$ (78) $\displaystyle T_{ij}$ $\displaystyle=0,$ (79)

and its trace is

$\displaystyle T_{\phantom{ij}}$ $\displaystyle\simeq-\rho.$ (80)

Inserting these relations into Eqs. (10) and (11), we obtain the system of
equations

$\displaystyle\vec{\nabla}\cdot\left(\varphi\vec{\nabla}\Phi\right)$
$\displaystyle\simeq 8\pi G\rho\,,$ (81)
$\displaystyle\varphi\nabla^{2}(\Phi-4U)-\vec{\nabla}U\cdot\vec{\nabla}\varphi\qquad\qquad$
$\displaystyle+\nabla^{2}\varphi-\frac{1}{\varphi}\vec{\nabla}\varphi\cdot\vec{\nabla}\varphi$
$\displaystyle\simeq 0\,,$ (82)
$\displaystyle\vec{\nabla}(U-\Phi)\cdot\vec{\nabla}\varphi-\nabla^{2}\varphi+2\Xi\varphi$
$\displaystyle\simeq 8\pi G\rho\,.$ (83)

Equation (81) derives from the $00$-component: it coincides with the modified
Poisson equation of RG when we identify $\varphi/2$ with the permittivity
$\epsilon$. Equations (82) and (83) derive from the contracted $ij$-components
and the scalar field equation, respectively. Given a density distribution
$\rho$, this system of equations can be solved for the three unknowns, $\Phi$,
$U$, and $\varphi$.

###  A.1 The gravitational field of a spherical source immersed in a
background of constant density

For spherically-symmetric systems, the exact solution of Eq. (81) is

$\displaystyle\frac{d\Phi}{dr}$
$\displaystyle=\frac{2}{\varphi}\frac{G\int_{0}^{r}\rho\,4\pi{r^{\prime}}^{2}dr^{\prime}}{r^{2}}$
$\displaystyle=\frac{2}{\varphi}\frac{Gm(<r)}{r^{2}}\,,$ (84)

where $\rho(r)=\rho_{s}(r)+\rho_{\mathrm{bg}}$, with $\rho_{s}(r)$ the density
profile of the source in a homogeneous background with density
$\rho_{\mathrm{bg}}$, and $m(<r)$ the enclosed mass within radius $r$.
Newtonian gravity is recovered when $\varphi=2$ (i.e. $\epsilon=1$). For
spherically-symmetric systems, Eqs. (82) and (83) are less straightforward to
solve than Eq. (81). In the following, we make a few approximations to derive
the generic behaviour of the gravitational field $d\Phi/dr$ and of the scalar
field $\varphi$ in the two limits of small and large distances from the
spherical source.

####  A.1.1 The gravitational field at small distances from a spherical source

In the Newtonian regime where $\rho_{s}\gg\rho_{\mathrm{bg}}$, we may
approximate $\varphi\simeq 2-\varphi_{1}$, where $|\varphi_{1}|\ll 1$. By
ignoring all the terms of the order $\mathcal{O}(\varphi_{1}^{2})$,
$\mathcal{O}(\varphi_{1})\mathcal{O}(\Phi)$, and
$\mathcal{O}(\varphi_{1})\mathcal{O}(U)$, Eq. (83) yields

$\nabla^{2}\varphi_{1}=\frac{1}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{d\varphi_{1}}{dr}\right)\simeq
8\pi G\rho\,,$ (85)

where we also ignored the terms containing $\Xi$: as we show in Sect. IV, we
can identify $\Xi$ with the cosmological constant of the standard model and,
due to its small value, we can thus safely ignore these terms in the weak-
field limit.

Solving Eq. (85) yields the scalar field

$\varphi(r)\simeq 2\left[1-\int_{0}^{r}\frac{Gm(<r^{\prime})}{r^{\prime
2}}\,dr^{\prime}\right]\,.$ (86)

If the spherical source has finite size $R$ and finite mass $m$, Eq. (86)
yields, beyond $R$,

$\varphi(r>R)\simeq 2+\frac{2Gm}{r}-\frac{4\pi
G\rho_{\mathrm{bg}}r^{2}}{3}\,,$ (87)

and Eq. (84) reads

$\frac{d\Phi}{dr}\simeq\frac{2Gm}{\varphi
r^{2}}\left(1+\frac{4\pi\rho_{\mathrm{bg}}r^{3}}{3m}\right)\,.$ (88)

In the vacuum, $\rho_{\mathrm{bg}}=0$, the RG field reduces to the Newtonian
field for $\varphi=2$, in agreement with Eq. (81), and it is larger than the
Newtonian field for $\varphi<2$.

####  A.1.2 The gravitational field at large distances from a spherical source

If the mass density profile of the extended spherical source decreases with
increasing $r$, at sufficiently large distances from the source, the
background density $\rho_{\mathrm{bg}}$ approximates the mean density of the
Universe, assuming that the source is isolated. We can thus treat the
contribution of the source, $\rho_{s}$, as a small perturbation to the total
density field $\rho=\rho_{\mathrm{bg}}+\rho_{s}$, with
$\rho_{s}\ll\rho_{\mathrm{bg}}$. We can similarly assume that the scalar field
reaches a mean cosmic value $\varphi_{\mathrm{bg}}=\mathrm{const}$ and write
the scalar field as $\varphi=\varphi_{\mathrm{bg}}+\varphi_{1}$, with
$\varphi_{1}\ll\varphi_{\mathrm{bg}}$. In the limit $\rho_{s}\to 0$ and
$\varphi_{1}\to 0$, Eqs. (82) and (83) simplify to

$\displaystyle\nabla^{2}U$
$\displaystyle\simeq\frac{1}{4}\nabla^{2}\Phi\simeq\frac{2\pi
G\rho_{\mathrm{bg}}}{\varphi_{\mathrm{bg}}}\,,$ (89)
$\displaystyle\varphi_{\mathrm{bg}}$ $\displaystyle\simeq\frac{4\pi
G\rho_{\mathrm{bg}}}{\Xi}\,.$ (90)

We can derive how the scalar field approaches the cosmic constant value
$\varphi_{\mathrm{bg}}$ by assuming that Eq. (89) already holds at distances
where $\varphi$ is not yet asymptotically constant, i.e. $\nabla^{2}U\simeq
2\pi G\rho/\varphi$. By adding up Eqs. (82) and (83), and using Eq. (81), we
obtain

$\displaystyle\left(\vec{\nabla}\ln\varphi\right)^{2}+2\vec{\nabla}\Phi\cdot\vec{\nabla}\ln\varphi-2\Xi+\frac{8\pi
G\rho}{\varphi}$ $\displaystyle\simeq 0\,.$ (91)

In spherical symmetry, we can recast this equation as

$\frac{d\ln\varphi}{dr}\simeq\frac{d\Phi}{dr}\left\\{-1\pm\left[1+\left(\frac{d\Phi}{dr}\right)^{-2}\left(2\Xi-\frac{8\pi
G\rho}{\varphi}\right)\right]^{1/2}\right\\}\,,$ (92)

which yields an implicit expression for the gravitational field at large
distances from the source. Hereafter, we consider the above equation with the
minus sign, so that $d\varphi/dr<0$.

We now explore the two limits of the large and small gravitational field
$d\Phi/dr$ at large distances from the source. In the limit
$d\Phi/dr\gg(2\Xi-{8\pi G\rho}/{\varphi})^{1/2}$, Eq. (92) reduces to
$d\varphi/dr\propto-2\varphi\,d\Phi/dr$, a dependence broadly approaching the
result derived in the previous section for small distances from the source
(Eqs. 84 and 86).

In the limit $d\Phi/dr\ll(2\Xi-{8\pi G\rho}/{\varphi})^{1/2}$, Eq. (92)
reduces to

$\frac{d\ln\varphi}{dr}\simeq-(2\Xi)^{1/2}\left(1-\frac{4\pi
G\rho}{\Xi\varphi}\right)^{1/2}$ (93)

which shows that $\varphi=\mathrm{const}$ when $\varphi_{\mathrm{bg}}=4\pi
G\rho/\Xi$, and thus $\rho=\rho_{\mathrm{bg}}$.

The relevant result of this analysis is that the transition between the two
regimes of large and small gravitational field, $d\Phi/dr$, takes place at the
acceleration scale

$\frac{d\Phi}{dr}\sim a_{\Xi}\equiv\left(2\Xi-\frac{8\pi
G\rho}{\varphi}\right)^{1/2}\,.$ (94)

We discuss this result in Sect. III.3.1.

At distances where $\varphi$ is not yet constant, namely
$\varphi(r)=\varphi_{\mathrm{bg}}+\varphi_{1}(r)$ and
$\rho(r)=\rho_{\mathrm{bg}}+\rho_{s}(r)$, we can write Eq. (93) as

$\frac{d\ln\varphi}{dr}\simeq-(2\Xi)^{1/2}\left(\frac{\varphi_{1}}{\varphi_{\mathrm{bg}}}-\frac{\rho_{s}}{\rho_{\mathrm{bg}}}\right)^{1/2}\,.$
(95)

Solving Eq. (95) for $\varphi$ requires an assumption on the form of
$\rho_{s}(r)$. For example, we can assume that the density profile
$\rho_{s}(r)$ of the source drops exponentially, i.e.
$\rho_{s}(r)=\rho_{\mathrm{edge}}e^{-r/r_{\mathrm{edge}}}$, where
$r_{\mathrm{edge}}$ is a characteristic scale of the source and
$\rho_{\mathrm{edge}}$ an appropriate normalization constant. In this case,
$\varphi_{1}(r)=\varphi_{\mathrm{edge}}e^{-r/r_{\mathrm{edge}}}$, with
$\varphi_{\mathrm{edge}}=4\pi G\rho_{\mathrm{edge}}/\Xi$, is a solution of Eq.
(95), if we ignore higher-order terms. We thus obtain the scalar field

$\varphi\simeq\frac{4\pi
G}{\Xi}\left(\rho_{\mathrm{bg}}+\rho_{\mathrm{edge}}e^{-r/r_{\mathrm{edge}}}\right)\,.$
(96)

##  Appendix B Derivation of the field equations for a static, spherically-
symmetric metric

Based on the metric of Eq. (27), the Christoffel symbols are777Throughout this
section we only show non-zero components.

$\displaystyle\Gamma^{0}_{r0}=\frac{\partial_{r}A}{A}\,,\quad\Gamma^{r}_{00}=\frac{A\partial_{r}A}{B^{2}}\,,\quad\Gamma^{r}_{rr}=\frac{\partial_{r}B}{B}\,,$
(97)
$\displaystyle\Gamma^{r}_{\theta\theta}=-\frac{r}{B^{2}}\,,\quad\Gamma^{r}_{\phi\phi}=-\frac{r\sin^{2}\theta}{B^{2}}\,,$
(98) $\displaystyle\Gamma^{\theta}_{\theta r}=\Gamma^{\phi}_{\phi
r}=\frac{1}{r}\,,\quad\Gamma^{\theta}_{\phi\phi}=-\sin\theta\cos\theta\,,$
(99)
$\displaystyle\Gamma^{\phi}_{\phi\theta}=\frac{\cos\theta}{\sin\theta}\,,$
(100)

from which we can calculate the components of the Ricci tensor

$\displaystyle R_{00}$
$\displaystyle=-\frac{A}{B}\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}A\right)\,,$
(101) $\displaystyle R_{rr}$
$\displaystyle=\frac{B}{A}\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}A\right)+2r\vec{\nabla}\cdot\left[\ln{(AB)}\vec{\nabla}\frac{1}{r}\right]\,,$
(102) $\displaystyle R_{\theta\theta}$
$\displaystyle=-\frac{r^{3}}{B^{2}}\vec{\nabla}\cdot\left(\ln{\frac{A}{B}}\vec{\nabla}\frac{1}{r}\right)+\frac{1}{B^{2}}-1\,,$
(103) $\displaystyle R_{\phi\phi}$ $\displaystyle=\sin^{2}\theta
R_{\theta\theta}\,.$ (104)

The three terms $\nabla_{\mu}\nabla_{\nu}\varphi$,
$\nabla^{\alpha}\nabla_{\alpha}\varphi$, and
$\nabla_{\mu}\varphi\nabla_{\nu}\varphi$ are

$\displaystyle\nabla_{0}\nabla_{0}\varphi$
$\displaystyle=-\frac{A}{B^{2}}\vec{\nabla}A\cdot\vec{\nabla}\varphi\,,$ (105)
$\displaystyle\nabla_{r}\nabla_{r}\varphi$
$\displaystyle=B\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}\varphi\right)+2r\vec{\nabla}\cdot\left(\varphi\vec{\nabla}\frac{1}{r}\right)\,,$
(106) $\displaystyle\nabla_{\theta}\nabla_{\theta}\varphi$
$\displaystyle=-\frac{r^{3}}{B^{2}}\vec{\nabla}\cdot\left(\varphi\vec{\nabla}\frac{1}{r}\right)\,,$
(107) $\displaystyle\nabla_{\phi}\nabla_{\phi}\varphi$
$\displaystyle=\sin^{2}\theta\nabla_{\theta}\nabla_{\theta}\varphi\,,$ (108)
$\displaystyle\nabla^{\alpha}\nabla_{\alpha}\varphi$
$\displaystyle=\frac{1}{AB}\vec{\nabla}\cdot\left(\frac{A}{B}\vec{\nabla}\varphi\right)\,,$
(109) $\displaystyle\nabla_{r}\varphi\nabla_{r}\varphi$
$\displaystyle=\vec{\nabla}\varphi\cdot\vec{\nabla}\varphi\,.$ (110)

For a perfect fluid, $T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu}$, we have

$\displaystyle T_{00}$ $\displaystyle=A^{2}\rho\,,$ (111) $\displaystyle
T_{rr}$ $\displaystyle=B^{2}p\,,$ (112) $\displaystyle T_{\theta\theta}$
$\displaystyle=r^{2}p\,,$ (113) $\displaystyle T_{\phi\phi}$
$\displaystyle=r^{2}\sin^{2}\theta p\,,$ (114) $\displaystyle T$
$\displaystyle=-\rho+3p\,.$ (115)

The field equations thus reduce to the following four equations

$\displaystyle\frac{1}{AB}\vec{\nabla}\cdot\left(\frac{\varphi}{2}\frac{1}{B}\vec{\nabla}A\right)=4\pi
G\rho\,,$ (116)
$\displaystyle\frac{\varphi}{2}\frac{1}{AB}\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}A\right)+\frac{\varphi}{2}\frac{2r}{B^{2}}\vec{\nabla}\cdot\left(\ln{(AB)}\vec{\nabla}\frac{1}{r}\right)$
$\displaystyle\qquad+\frac{1}{B}\vec{\nabla}\cdot\left(\frac{1}{B}\vec{\nabla}\frac{\varphi}{2}\right)+\frac{2r}{B^{2}}\vec{\nabla}\cdot\left(\frac{\varphi}{2}\vec{\nabla}\frac{1}{r}\right)$
$\displaystyle\qquad\qquad-\frac{2}{\varphi
B^{2}}\vec{\nabla}\frac{\varphi}{2}\cdot\vec{\nabla}\frac{\varphi}{2}=-4\pi
Gp\,,$ (117)
$\displaystyle\frac{\varphi}{2}\left[-\frac{r}{B^{2}}\vec{\nabla}\cdot\left(\ln{\frac{A}{B}}\vec{\nabla}\frac{1}{r}\right)+\frac{1}{r^{2}B^{2}}-\frac{1}{r^{2}}\right]$
$\displaystyle\qquad\qquad-\frac{r}{B^{2}}\vec{\nabla}\cdot\left(\frac{\varphi}{2}\vec{\nabla}\frac{1}{r}\right)=-4\pi
Gp\,,$ (118)
$\displaystyle\frac{1}{2AB}\vec{\nabla}\cdot\left(\frac{A}{B}\vec{\nabla}\varphi\right)-\Xi\varphi=4\pi
G(-\rho+3p)\,.$ (119)

##  Appendix C Derivation of the modified Friedmann equations and of the
scalar field equation

Based on the FLRW metric, Eq. (33), the Christoffel symbols are888Throughout
this section we only show non-zero components.

$\displaystyle\Gamma^{0}_{rr}=\frac{a\dot{a}}{1-kr^{2}}\,,\quad\Gamma^{0}_{\theta\theta}=a\dot{a}r^{2}\,,$
(120)
$\displaystyle\Gamma^{0}_{\phi\phi}=a\dot{a}r^{2}\sin^{2}\theta\,,\quad\Gamma^{r}_{0r}=\frac{\dot{a}}{a}\,,\quad\Gamma^{\theta}_{0\theta}=\Gamma^{\phi}_{0\phi}=\frac{\dot{a}}{a}\,,$
(121)
$\displaystyle\Gamma^{r}_{rr}=\frac{kr}{1-kr^{2}}\,,\quad\Gamma^{r}_{\theta\theta}=-r(1-kr^{2})\,,$
(122)
$\displaystyle\Gamma^{r}_{\phi\phi}=-r(1-kr^{2})\sin^{2}\theta\,,\quad\Gamma^{\theta}_{r\theta}=\Gamma^{\phi}_{r\phi}=\frac{1}{r}\,,$
(123)
$\displaystyle\Gamma^{\theta}_{\phi\phi}=-\sin\theta\cos\theta\,,\quad\Gamma^{\phi}_{\theta\phi}=\frac{\cos\theta}{\sin\theta}\,,$
(124)

The components of the Ricci tensor are

$\displaystyle R_{00}$ $\displaystyle=3\frac{\ddot{a}}{a}\,,$ (125)
$\displaystyle R_{rr}$
$\displaystyle=-\frac{1}{1-kr^{2}}\left(a\ddot{a}+2\dot{a}^{2}+2k\right)\,,$
(126) $\displaystyle R_{\theta\theta}$
$\displaystyle=-r^{2}\left(a\ddot{a}+2\dot{a}^{2}+2k\right)\,,$ (127)
$\displaystyle R_{\phi\phi}$
$\displaystyle=-r^{2}\sin^{2}\theta\left(a\ddot{a}+2\dot{a}^{2}+2k\right)\,,$
(128)

and the three terms $\nabla_{\mu}\nabla_{\nu}\varphi$,
$\nabla^{\alpha}\nabla_{\alpha}\varphi$, and
$\nabla_{\mu}\varphi\nabla_{\nu}\varphi$ are

$\displaystyle\nabla_{0}\nabla_{0}\varphi=\ddot{\varphi}\,,\quad\nabla_{r}\nabla_{r}\varphi=-\frac{1}{1-kr^{2}}a\dot{a}\dot{\varphi}\,,$
(129)
$\displaystyle\nabla_{\theta}\nabla_{\theta}\varphi=-r^{2}a\dot{a}\dot{\varphi}\,,\quad\nabla_{\phi}\nabla_{\phi}\varphi=-r^{2}\sin^{2}\theta
a\dot{a}\dot{\varphi}\,,$ (130)
$\displaystyle\nabla^{\alpha}\nabla_{\alpha}\varphi=-\ddot{\varphi}-3\frac{\dot{a}}{a}\dot{\varphi}\,,\quad\nabla_{0}\varphi\nabla_{0}\varphi=\dot{\varphi}^{2}.$
(131)

From the stress-energy tensor of a perfect fluid,
$T_{\mu\nu}=\left(\rho+p\right)u_{\mu}u_{\nu}+pg_{\mu\nu}$, we get

$\displaystyle T_{00}=\rho\,,\quad T_{rr}=\frac{a}{1-kr^{2}}p\,,\quad
T_{\theta\theta}=r^{2}a^{2}p\,,$ (132) $\displaystyle
T_{\phi\phi}=r^{2}\sin^{2}\theta a^{2}p\,,\quad T=-(\rho-3p)\,.$ (133)

By combining the above results, the time-time component of the modified
Einstein field equations, Eq. (10), is

$\varphi
R_{00}+\nabla_{0}\nabla_{0}\varphi-\frac{1}{\varphi}\nabla_{0}\varphi\nabla_{0}\varphi=-8\pi
GT_{00}\,,$

which gives the first modified Friedmann equation

$\frac{\ddot{a}}{a}-\frac{1}{3}\left(\frac{\dot{\varphi}^{2}}{\varphi^{2}}-\frac{\ddot{\varphi}}{\varphi}\right)=-\frac{8\pi
G}{3\varphi}\rho\,.$ (134)

From the $rr$-component of Eq. (10),

$\varphi
R_{rr}+\nabla_{r}\nabla_{r}\varphi-\frac{1}{\varphi}\nabla_{r}\varphi\nabla_{r}\varphi=-8\pi
GT_{rr}\,,$ (135)

we obtain the second modified Friedmann equation

$\frac{\ddot{a}}{a}+2\frac{\dot{a}^{2}}{a^{2}}+\frac{2k}{a^{2}}+\frac{\dot{a}\dot{\varphi}}{a\varphi}=\frac{8\pi
G}{\varphi}p\,.$ (136)

The $\theta\theta$\- and $\phi\phi$-components reduce to the same expression.

The last equation is the scalar field equation

$\Box\varphi-2\Xi\varphi=8\pi GT\,,$ (137)

which becomes

$\frac{1}{3}\frac{\ddot{\varphi}}{\varphi}+\frac{\dot{a}\dot{\varphi}}{a\varphi}+\frac{2\Xi}{3}=\frac{8\pi
G}{3\varphi}\left(\rho-3p\right)\,.$ (138)

##  Appendix D Derivation of the time derivative of the Hubble parameter [Eq.
(56)]

Adding up Eqs. (49) and (51) yields

$q=\Omega-2\Omega_{\Xi}-\zeta-\frac{1}{3}\zeta^{2}\,.$ (139)

Using the definitions of $q$ and $\zeta$, Eqs. (44) and (48), together with
Eq. (54), the above equation reduces to

$\frac{\ddot{a}^{2}a^{2}}{2\dot{a}^{4}}-\frac{\ddot{a}a}{\dot{a}^{2}}=\frac{3}{2}\Omega-3\Omega_{\Xi}+1\,.$
(140)

With Eqs. (54) and (38) and the expression
$\Omega(t_{0})\equiv\Omega_{0}=16\pi G\rho_{0}/(3H_{0}^{2}\varphi_{0})$, we
can recast Eq. (42) as

$\Omega=\frac{\Omega_{0}H_{0}}{H}\,,$ (141)

and similarly (Eq. 43),

$\Omega_{\Xi}=\frac{\Omega_{\Xi 0}H_{0}^{2}}{H^{2}}\,.$ (142)

Equation (140) thus becomes

$\frac{\ddot{a}^{2}}{\dot{a}^{2}}-\frac{2\ddot{a}}{a}-\frac{2\dot{a}^{2}}{a^{2}}=3\Omega_{0}H_{0}\frac{\dot{a}}{a}-6\Omega_{\Xi
0}H_{0}^{2}\,,$ (143)

and, by adding $3\dot{a}^{2}/a^{2}$ on both sides, we obtain

$\left(\frac{\ddot{a}}{\dot{a}}-\frac{\dot{a}}{a}\right)^{2}=3\frac{\dot{a}^{2}}{a^{2}}+3\Omega_{0}H_{0}\frac{\dot{a}}{a}-6\Omega_{\Xi
0}H_{0}^{2}\,.$ (144)

The right-hand side can be rewritten as

$\left(\frac{d}{dt}\ln\dot{a}-\frac{d}{dt}\ln
a\right)^{2}=\left(\frac{d}{dt}\ln\frac{\dot{a}}{a}\right)^{2}\,,$ (145)

which leads to the diffential equation for $H$

$\left(\frac{dH}{dt}\right)^{2}=3\left(H^{4}+\Omega_{0}H_{0}H^{3}-2\Omega_{\Xi
0}H_{0}^{2}H^{2}\right)\,.$ (146)

##  Appendix E Analytic derivation of the luminosity distance

The general definition of the luminosity distance is

$D_{L}=(1+z)\int_{t}^{t_{0}}\frac{dt^{\prime}}{a(t^{\prime})}.$ (147)

We will compute its exact form using the solution of the scale factor in Eq.
(59), that we rewrite here as

$\begin{split}a^{\sqrt{3}}C_{2}^{-\sqrt{3}}=\left|1-\frac{2\sqrt{8\Omega_{\Xi
0}}}{\Omega_{0}\tan\left[\sqrt{\frac{3}{2}\Omega_{\Xi
0}}\left(H_{0}t+\mathcal{C}_{1}\right)\right]\mp\sqrt{\Omega_{0}^{2}+8\Omega_{\Xi
0}}+\sqrt{8\Omega_{\Xi 0}}}\right|.\end{split}$ (148)

Inverting Eq. (148) yields

$\begin{split}&\sqrt{\frac{3}{2}\Omega_{\Xi
0}}\left(H_{0}t+C_{1}\right)\equiv\arctan\left[A+\frac{B}{1-a^{\sqrt{3}}C_{2}^{-\sqrt{3}}}\right],\end{split}$
(149)

where we defined

$A\equiv\pm\frac{\sqrt{\Omega_{0}^{2}+8\Omega_{\Xi
0}}}{\Omega_{0}}-\frac{\sqrt{8\Omega_{\Xi 0}}}{\Omega_{0}},\quad
B\equiv\frac{2\sqrt{8\Omega_{\Xi 0}}}{\Omega_{0}}.$ (150)

Differentiating Eq. (149) yields

$\begin{split}&dt=\frac{8}{H_{0}\Omega_{0}}\frac{C_{2}^{-\sqrt{3}}a^{\sqrt{3}-1}}{\left(1-a^{\sqrt{3}}C_{2}^{-\sqrt{3}}\right)^{2}+\left(A+B-Aa^{\sqrt{3}}C_{2}^{-\sqrt{3}}\right)^{2}}da.\end{split}$
(151)

The luminosity distance in Eq. (147) then becomes

$\begin{split}D_{L}=\frac{8(1+z)}{H_{0}\Omega_{0}}\int_{0}^{z}\frac{C_{2}^{-\sqrt{3}}(1+z)^{-\sqrt{3}}}{\left[1-\left(1+z\right)^{-\sqrt{3}}C_{2}^{-\sqrt{3}}\right]^{2}+\left[A+B-A\left(1+z\right)^{-\sqrt{3}}C_{2}^{-\sqrt{3}}\right]^{2}}dz.\end{split}$
(152)

If we change variable by defining

$\begin{split}&y\equiv\left(1+z\right)^{-\sqrt{3}}C_{2}^{-\sqrt{3}}\,,\end{split}$
(153)

the luminosity distance becomes

$D_{L}=\frac{8\left(1+z\right)}{\sqrt{3}C_{2}H_{0}\Omega_{0}}\int_{y}^{\mathcal{C}_{2}^{-\sqrt{3}}}\frac{y^{-1/\sqrt{3}}dy}{(1-y)^{2}+(A+B-Ay)^{2}}\,,$
(154)

with

$\begin{split}&\int\frac{y^{-1/\sqrt{3}}dy}{(1-y)^{2}+(A+B-Ay)^{2}}=\frac{i\sqrt{3}}{2B}y^{-1/\sqrt{3}}\biggl{\\{}\frac{(A-i)y}{(A-i)(y-1)-B}\
_{2}F_{1}\left(1,1;1+\frac{1}{\sqrt{3}};\frac{A+B-i}{(A-i)(1-y)+B}\right)\\\
&-\frac{(A+i)y}{(A+i)(y-1)-B}\
_{2}F_{1}\left(1,1;1+\frac{1}{\sqrt{3}};\frac{A+B+i}{(A+i)(1-y)+B}\right)\biggr{\\}}.\end{split}$
(155)

Equation (154) corresponds to two equations, according to the sign of the
parameter $A$ defined in Eq. (150). The two curves corresponding to Eq. (154)
appear in Fig. 1.