Covariant formulation of refracted gravity

We propose a covariant formulation of refracted gravity (RG), a classical theory of gravity based on the introduction of gravitational permittivity (GP), a monotonic function of the local mass density, in the standard Poisson equation. GP mimics dark matter (DM) phenomenology. The covariant formulation of RG (CRG) that we propose belongs to the class of scalar-tensor theories, where the scalar field $\varphi$ has a self-interaction potential $V(\varphi)=-\Xi\varphi$, with $\Xi$ a normalization constant. We show that $\varphi$ is twice the GP 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_{bg}$, $\rho_{bg}$ being 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 is comparable to the acceleration $a_0$ originally introduced in MOND. From CRG, we also derived the modified Friedmann equations for an expanding, homogeneous, and isotropic universe. We find that the same scalar field that mimics DM also drives the accelerated expansion of the Universe. From the stress-energy tensor of $\varphi$, we derived the equation of state of a redshift-dependent effective dark energy (DE) $w_{DE}=p_{DE}/\rho_{DE}$. Current observational constraints on $w_{DE}$ and distance modulus data of SNIa suggest that $\Xi$ has a comparable value to the cosmological constant $\Lambda$ in the standard model. CRG, therefore, suggests a natural explanation of the known relation $a_0\sim \Lambda^{1/2}$ and appears to describe both the dynamics of cosmic structure and the expanding Universe with a single scalar field, highlighting a possible deep connection between phenomena currently attributed to DM and DE separately.


Introduction
The current standard Λ cold dark matter (ΛCDM) cosmological model assumes that gravitational interactions are ruled by general relativity (GR); the model relies on the existence of collisionless non-baryonic CDM, and a positive cosmological constant Λ (Ostriker & Steinhardt 1995). Non-baryonic dark matter is required to account for the abundance of light elements (Cyburt et al. 2016), the amplitude of the power spectrum of the temperature anisotropies in the cosmic microwave background (CMB) (Aghanim et al. 2020b), and the dynamics of cosmic structures (Clowe et al. 2006;Dodelson et al. 2001;Del Popolo 2014;Akrami et al. 2020;Aghanim et al. 2020a). The cosmological density parameter Ω Λ0 = Λ/3H 2 0 associated with Λ, with H 0 being the Hubble constant, accounts for the negative value of the deceleration parameter q 0 measured with the Hubble diagram of type Ia supernovae (SNeIa) (Riess et al. 1998;Perlmutter et al. 1999). The curvature of the Universe, Ω k = Ω Λ0 + Ω 0 − 1, measured from the CMB power spectrum, suggests a null curvature (Aghanim et al. 2020c).
Although the Λ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 (Riess et al. 2016) is more than 4σ larger than H 0 inferred from the measurements of the CMB (Verde et al. 2019). The observed abundances of the light elements are consistent with the Big Bang nucleosynthesis scenario, except 7 Li, whose abundance is a factor of ∼ 3 smaller (Mathews et al. 2020). Numerous features of the CMB temperature anistropies are present on large scales (Ade et al. 2016). The probability of some of these features to appear individually is 0.1%; the combined probability of the uncorrelated features is 0.001% and might represent a serious challenge to the ΛCDM model (Schwarz et al. 2016;Luongo et al. 2022). The lensing amplitude in the CMB power spectra is enhanced compared to ΛCDM expectations (Aghanim et al. 2020b) and would suggest a positive rather than a null curvature of the Universe (Handley 2021;Di Valentino et al. 2019). A slight tension also appears for the normalization of the power spectrum σ 8 (e.g. Raveri 2016): the value inferred from the CMB measurements (Aghanim et al. 2020c) is 2σ larger than the value derived from the tomographic weak gravitational lensing analysis of the Kilo-Degree Survey (KiDS) imaging data (Hildebrandt et al. 2017).
In addition, the cosmological constant poses a fine-tuning problem that is theoretically challenging. If we associate Λ with the ground state energy level of the vacuum in quantum field theory, its measured value, Λ ∼ 10 −12 eV 4 , appears to be ∼ 120 orders of magnitude smaller than expected (Weinberg 1989;Padilla 2015). If we associate Λ with the energy scale up to which the standard model of particle physics has been tested, ∼ 1 TeV, the discrepancy reduces to ∼ 60 orders of magnitude  (Joyce et al. 2015), but it remains severe. The most popular solution to the Λ problem is to suppress Λ 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 & Ratra 2003;Copeland et al. 2006;Bamba et al. 2012;Joyce et al. 2015;Amendola et al. 2018, for extensive reviews).
Collisionless CDM poses additional problems on small scales: the core-cusp problem in dwarf and disk galaxies, the missing satellite problem, the too-big-to-fail problem, and the plane of satellite galaxies in the Milky Way and nearby large galaxies (Salucci 2003;Ferrero et al. 2012;Boylan-Kolchin et al. 2012;Garrison-Kimmel et al. 2014;Del Popolo & Le Delliou 2017;Kroupa 2012;de Martino et al. 2020). In addition, some relations, such as the radial acceleration relation or the baryonic Tully-Fisher relation in disk galaxies, would require finely tuned interactions between CDM and baryonic matter (Desmond & Wechsler 2015;Di Cintio & Lelli 2016;Desmond 2017).
Some of these small-scale tensions may originate either by an inaccurate treatment of the dynamics of CDM and baryonic matter or by the inappropriate properties adopted for the dark matter model. For example, the core-cusp problem emerges when we attempt to interpret the observed kinematics of stars in galaxies by assuming that the galaxies are embedded within CDM halos with a Navarro-Frenk-White (NFW) density profile, as predicted by N-body simulations (Navarro et al. 1997). By dropping this constraint on the dark matter distribution and adopting a dark matter density profile with a flat central core, we can properly describe the stellar kinematics of spirals with a universal rotation curve (Persic et al. 1996;Salucci et al. 2007;Salucci & De Laurentis 2012;Salucci 2018). Indeed, an exponential disk and a dark matter halo described by the Burkert profile with a core (Burkert 1995) excellently describe the rotation curves of five spirals (Gentile et al. 2004), and suggest that spirals and dwarf galaxies lie on the same scaling relation between the core density and the core radius of the dark matter halo (Salucci & Burkert 2000). A number of different effects from the dynamics of CDM or from baryonic physics, including stellar feedback and star formation efficiency, are advocated to generate a central core in the dark matter distribution (de Martino et al. 2020). For example, tidal effects by a massive hosting galaxy might induce dark matter density profiles shallower than the NFW profile in the central regions of satellite halos and might also alleviate the too-big-to-fail problem (Tomozeiu et al. 2016).
Dropping the hypothesis that dark matter is collisionless and cold might solve some, albeit not all, of the tensions of CDM on the scale of galaxies (Salucci 2019;Salucci & di Paolo 2021). Indeed, weakly interacting massive particles (WIMPs) are the most plausible candidate to make up the collisionless CDM (Bertone et al. 2005). However, attempts to, directly or indirectly, detect these particles have not yet produced unquestionable results (Tanabashi et al. 2018). Alternative dark matter models include warm dark matter, self-interacting dark matter, QCD axions, and fuzzy dark matter (de Martino et al. 2020). Some of these models have particle counterparts, such as sterile neutrinos or ultra-light bosons. The detection of these particles requires direct or indirect experiments different from those conceived for detecting WIMPs (Buckley & Peter 2018). Dark matter particles with distinct peculiar features, such as superfluidity (Berezhiani & Khoury 2015) or gravitational polarization (Blanchet 2007;Blanchet & Heisenberg 2017), can also partly reproduce the phenomenology of galaxies.
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 MOdified Newtonian Dynamics (MOND) was originally motivated only by the observations of the flat rotation curves of disk galaxies (Milgrom 1983a;Begeman et al. 1991;Begum & Chengalur 2004), and many of the current observational challenges on the scale of galaxies were actually predicted by MOND (Milgrom 1983b,c;Sanders & McGaugh 2002). The universal rotation curves of spirals (Persic et al. 1996) and the universality of the galactic surface density within the radius of the core of a Burkert dark matter profile also appear to be consistent with MOND (Gentile 2008;Gentile et al. 2009), although the debate on the dynamics of dwarf disk galaxies remains vibrant (Corbelli & Salucci 2007;Sanchez-Salcedo et al. 2013;Di Paolo et al. 2019;Banik et al. 2020;Salucci & di Paolo 2021). Furthermore, the MOND formulation is purely phenomenological and its extension to a covariant formulation has proven to be challenging (Famaey & McGaugh 2012;Milgrom 2015;Skordis & Złośnik 2019;Złośnik & Skordis 2017;Skordis & Zlosnik 2021).
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. (2019) 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 1991;Padmanabhan 2015;Verlinde 2017), 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 et al. 2018a,b;Tuveri & Cadoni 2019;Cadoni et al. 2020).
Attempts to describe both the accelerated expansion of the Universe and the dynamics of cosmic structures without dark matter and dark energy by building a modified theory of gravity are also numerous (see, e.g., Clifton et al. 2012;Nojiri et al. 2017). According to Lovelock's theorem (Lovelock 1971(Lovelock , 1972, 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 & Kazanas 1994). Conformal gravity does not present ghost-like instabilities, that might be common in theories with high-order derivatives (Bender & Mannheim 2008), and is renormalizable (Mannheim 2012). Unfortunately, although conformal gravity successfully reproduces the accelerated expansion of the Universe , the expected abundance of primordial deuterium is orders of magnitudes smaller than observed (Elizondo & Yepes 1994). Moreover, conformal gravity is unable to reproduce the dynamics of galaxy clusters (Horne 2006;Diaferio & Ostorero 2009), and appears to require a fine-tuning condition to describe the phenomenology of gravitational lensing and the dynamics of disk galaxies (Campigotto et al. 2019).
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 2018). 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 et al. 2010). These models can be a viable alternative to ΛCDM (Camera et al. 2011(Camera et al. , 2019 if the effective sound speed is small enough that the scalar field can cluster (Bertacca et al. 2008;Camera et al. 2009).
Here, we contribute to the quest for a unified model of the two dark sectors by proposing a novel scalar-tensor theory (Fujii & Maeda 2007;Quiros 2019) 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), which is a new phenomenological modified theory of gravity proposed by Matsakos & Diaferio (2016). Refracted gravity 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. (2020) showed that RG properly describes the rotation curves and the vertical velocity dispersion profiles of 30 disk galaxies in the DiskMass Survey (Bershady et al. 2010), and the dynamics of stars and globular clusters in the outer regions of three elliptical galaxies of type E0 (Cesare et al. 2022). Here, we provide a covariant formulation of this non-relativistic formulation of RG.
Section 2 reviews the relevant features of RG. In Sect. 3, 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. 4, 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 derive the equation of state of a redshift-dependent effective dark energy. We conclude in Sect. 5.

Non-relativistic refracted gravity
The phenomenological RG is based on the modified Poisson equation (hereafter RG equation) (Matsakos & Diaferio 2016) where Φ is the gravitational potential, G the gravitational constant, and ρ the density of ordinary matter. The function 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, that is (ρ). However, ρ was only chosen because is the simplest scalar field characterizing the matter distribution in the weak-field regime; 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, must be constant and equal to 1 in order to recover the standard Newtonian Poisson equation, ∇ 2 Φ = 4πGρ. Therefore, we can adopt the form of the permittivity (ρ) 1 for ρ ρ thr v for ρ ρ thr , where ρ thr is the density threshold that sets the transition between the Newtonian regime and the modified gravity regime; 1 0 < v < 1 is the permittivity of the vacuum. The presence of 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; and (2) it bends the gravitational field lines in regions where ∇Φ and ∇ are not parallel. The former effect is trivially seen in the simple case = const < 1, when ρ < ρ thr : Eq. (1) becomes ∇ 2 Φ = 4πGρ/ , and shows that the resulting field is equivalent to a Newtonian field originating from a larger effective mass density ρ/ , or a larger gravitational constant G/ . When ∇Φ and ∇ are not parallel, the field lines are refracted; in other words, the field lines change their direction when they cross the iso-surfaces of at an angle different from π/2. This effect can also generate non-Newtonian phenomena in regions where the density is larger than ρ thr , because the redirection of the field lines has nonlocal 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 ρ > ρ thr . The redirection of the field lines is also expected to be responsible for the mass discrepancy in galaxy clusters (Matsakos & Diaferio 2016). In these two studies, (ρ) is assumed to be a smooth step function, with the density threshold ρ thr in the range 10 −27 -10 −24 g cm −3 .

Refracted gravity as a scalar-tensor theory
We present the fundamental equations of CRG in Sect. 3.1, and the relation of CRG with the Horndeski theory and the CRG screening mechanisms in Sect. 3.2 and Sect. 3.3, respectively. In Sect. 3.4, we derive the weak field limit of CRG.

Fundamental equations of CRG
The family of scalar-tensor (TeS) theories, with the scalar field ϕ non-minimally coupled to the rank-2 tensor field g µν , can be derived from the action 2 where −g is the determinant of g µν , R = g µν R µν is the Ricci scalar, R µν is the Ricci tensor, g µν is the inverse metric, ∇ µ is the covariant derivative, ∇ α ϕ∇ α ϕ ≡ g αβ ∇ α ϕ∇ β ϕ, and L m is the matter Lagrangian density, with ψ m being the matter fields (Faraoni 2004). The potential V(ϕ) and the general differentiable function of the scalar field W(ϕ) parametrize the family of TeS theories. When ϕ → 1, the TeS action reduces to the standard Einstein-Hilbert action with a cosmological constant equal to −V(1). Hereafter, for the sake of simplicity, we consider the ϕdependence of W and V implicit.
(3) with respect to the scalar field yields the equation for ϕ We define CRG by setting We note that the original Brans-Dicke theory -with a constant W of order unity and a zero potential -is ruled out by post-Newtonian expansions and solar system experimental tests, which give the constraint W 40 000 (Faraoni 2004;Clifton et al. 2012), and by recent results from CMB observations (e.g. Li et al. 2013;Avilez & Skordis 2014). However, these constraints seem to be weaker on large cosmological scales and can be avoided by adding a self-interaction potential (Hrycyna et al. 2014;Quiros 2019), which we define as with Ξ being a constant parameter.
With the definitions of Eqs. (6) and (7), the modified field equations and the equation for ϕ, in the JF, become By using Eq. (9) to simplify Eq. (8), and by using Eq. (8), contracted with g µν , to simplify the resulting Eq. (9), we obtain the following CRG equations: We derive the stress-energy tensor for the scalar field by recasting Eq. (8) as 3 We adopt the sign conventions of Weinberg (1972), with the Ricci tensor given by , and the standard Einstein equations G µν = −8πGT µν . 4 We choose to work in the JF rather than in the Einstein frame (EF) because, in the JF, the scalar field is interpreted as a modification to the gravitational field (the left-hand side of the field equations) rather than a modification to the stress-energy tensor as in the EF (the right-hand side of the field equations). In the weak-field limit, this modification is consistent with the Poisson equation of the non-relativistic formulation of RG (Eq. 1), which contains modifications to the field, and not to the source term. As we show in Sect. 3.4, this formulation naturally leads to the identification of the scalar field with the permittivity.
where T M µν refers to the matter/energy stress-energy tensor, and T ϕ µν refers to the scalar-field stress-energy tensor, whose explicit form is 5 We use the stress-energy tensor of ϕ to derive an effective dark energy in Sect. 4.3.

Relation to Horndeski theories
Scalar-tensor theories belong to the wider family of Horndeski models (for reviews on the topic, see, e.g., Joyce et al. 2015;Kobayashi 2019), and, consequently, CRG can also be identified as a special case of Horndeski theories.

Screening mechanism of CRG
General relativity 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 et al. 2015).
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 advantage of the EF is that the field equations have a manifestly GR-like form 7 and computations can be performed more easily. However, the presence of the extra coupling between the scalar and matter fields alters rods and clocks, and makes the identification of physical observables harder.
The transition from the JF to the EF can be performed with the conformal transformation of the metric (Clifton et al. 2012), namely a length-scale transformation 8 g µν = e −2Γ g µν and g = e −4Γ √ g , where the tilde indicates quantities computed in the EF and Γ = − 1 2 ln ϕ. From the JF action of Eq. (3), the corresponding action in the EF is where Ψ is the scalar field related to the scalar field in the JF by Ψ = ∓ ln ϕ,R is the Ricci scalar, A(Ψ) ≡ ϕ −1/2 is the coupling between the scalar field and matter, and U(Ψ) is the potential (Clifton et al. 2012) The potential U(Ψ) = 2Ξe −Ψ is the runaway potential of the chameleon cosmology (Clifton et al. 2012;Joyce et al. 2015;Khoury & Weltman 2004). This potential can be studied by considering the non-relativistic equation of motion for the scalar field (see, e.g., Clifton et al. 2012, and references therein) where V eff is the effective potential. In particular, d 2 V eff /dΨ 2 can be interpreted as an effective mass of the scalar field, m Ψ , which reads as In regimes of large densities ρ e −2Ψ Ξ/2πG, the scalar field becomes increasingly more massive and hence the fifth force mediated by it has an increasingly shorter range; in other words, in this regime the presence of the scalar field is effectively screened. Conversely, in regimes of small densities ρ e −2Ψ Ξ/2πG, m Ψ gets smaller, the force mediated by Ψ has an increasingly longer range, and the field is free to propagate. 9 The density scale separating the two regimes thus depends on both Ψ and Ξ. This result scalar field acts as an additional gravitational source. Moreover, the extra coupling between the scalar and matter fields introduces both deviations from the geodesic motions of free-falling particles and a stressenergy tensor which is not covariantly conserved (Faraoni 2004;Clifton et al. 2012). 8 Since a direct measurement of an absolute scale is not possible, experiments are unable to distinguish between the EF and the JF frames. These frames represent two different realizations of the same theory, and physical observables must be equivalent in the two frames (see, e.g., Sotiriou et al. 2008;Postma & Volponi 2014). 9 A theorem guarantees that in TeS theories endowed with a potential U that satisfies the condition d 2 U(Ψ)/dΨ 2 > 0, like CRG (Eq. 17), black holes in vacuum are equivalent to GR black holes (Bekenstein 1995;Sotiriou & Faraoni 2012;Cruz et al. 2017). In CRG, black holes embedded in environment with density sufficiently small to make the screening mechanism ineffective, might in principle develop scalar hair. However, in extended theories of gravity, black holes or compact objects with scalar hair remain viable and their existence can be tested with gravitational wave observations (Sotiriou 2015;Berti et al. 2015;Cardoso et al. 2016;Brito et al. 2017;Barack et al. 2019;Maggiore et al. 2020). is consistent with the assumption, in the phenomenological RG, that the gravitational permittivity depends on the local density, and suggests a relation between and Ψ, or equivalently ϕ. In addition, this result indicates that the density scale depends on Ξ, which is a constant universal value independent of the local environment. In Sect. 4.3 and Appendix E, we find that Ξ plays the role of the cosmological constant Λ in the standard cosmological model. Therefore, we expect that the local value of the scalar field Ψ (or ϕ), rather than Ξ, plays the major role in setting the density threshold for the screening mechanism.
3.4. The weak-field limit of CRG As in standard GR, we take the static weak-field limit of the metric g µν , 10 which can be expanded around the Minkowski metric η µν as where |h| 1 and |∂ µ h| 1. We write the metric as where Φ and U are two potentials, and we ignore terms of order O(Φ 2 ) ∼ O(U 2 ). These expressions enable the calculation of the left-hand side of Eqs. (10) and (11) (10) and (11), we consider a static non-relativistic fluid: the only non-zero component of the energy-momentum tensor is T 00 ρ, and thus its trace is T −ρ. With these assumptions, the 00-component of the field equations reduces to ∇ · (ϕ∇Φ) 8πGρ .
We thus recover the RG equation if we identify the scalar field with twice the permittivity: ϕ = 2 . In the Newtonian regime, we have a constant scalar field, namely ∇ϕ = 0, and thus we recover the standard Poisson equation for ϕ = 2.
3.4.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 ρ bg . The source is described by a density profile ρ s (r) decreasing with r. We estimate the field in the limit ρ s (r) ρ bg , close to the source, and ρ s (r) ρ bg , at large distance from the source.
At small distances, we find with the scalar field At large distances from the source, the scalar field ϕ and the gravitational field dΦ/dr satisfy the implicit relation with ρ(r) = ρ s (r) + ρ bg . This relation sets the acceleration scale In regions where dΦ/dr a Ξ , the gravitational field dΦ/dr −(1/2)d ln ϕ/dr has a similar dependence of the field at small distances (Eqs. 25 and 26). In regions where dΦ/dr a Ξ , 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. 4, we find Ξ ∼ Λ. We thus find that a Ξ ∼ 10 −10 m s −2 , in the limit 2Ξ 8πGρ/ϕ, that occurs at large distances from the source. The existence of this acceleration scale appears in a number of observations on the scale of galaxies (see, e.g., de Martino et al. 2020;McGaugh 2020;Chae et al. 2020;McGaugh et al. 2016McGaugh et al. , 2018Merritt 2020). 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 et al. 2022). Nevertheless, some of the phenomenology predicted by MOND in the low-acceleration regime, set by an acceleration scale a 0 independent of the source, like the radial acceleration relation (McGaugh et al. 2016), appears to be inconsistent with the rotation curves of dwarf disk spirals and low-surface brightness galaxies (Di Paolo et al. 2019;Santos-Santos et al. 2020). These tensions might suggest that the acceleration scale could indeed depend on the source, as it happens for a Ξ in CRG.
The connection between a Ξ and Ξ, in the limit 2Ξ 8πGρ/ϕ, is similar to the connection between a 0 and Λ in MOND: a 0 ∼ Λ 1/2 . A number of different sensible arguments have been suggested for the interpretation of the latter relation (Milgrom 1989(Milgrom , 1999Famaey & McGaugh 2012;Milgrom 2020). In the CRG context, a Ξ ∼ (2Ξ) 1/2 emerges naturally.

The scalar field ϕ for a spherical source immersed in a constant background
According to the results of Appendix A.1, the scalar field ϕ is positive and broadly in the range [0, 2] consistently with the RG ansatz ∈ [0, 1] for the permittivity. In the limit ρ bg ρ s , we have whereas in the limit ρ bg ρ s , These limits are broadly in agreement with the smooth step function considered for = (ρ) in previous RG studies (Cesare et al. 2022;Cesare et al. 2020;Matsakos & Diaferio 2016). In those studies, the local mass density was found to be a good proxy of the transition between the Newtonian and the RG regimes.

The homogeneous and isotropic universe in CRG
We derive the basic equations of a homogeneous and isotropic universe in Sect. 4.1 and we solve these equations for a spatially flat universe in Sect. 4.2. We derive the equation of state of the effective dark energy in Sect. 4.3 and discuss the evolution of the scalar field in Sect. 4.4.

Basic equations
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 dr 2 1 − kr 2 + r 2 dθ 2 + r 2 sin 2 θ dφ 2 , where a(t) is the scale factor, k is the spatial curvature, and (t, r, θ, φ) are 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 B) and The equation for the scalar field, Eq. (11), reduces to (see Appendix B) In the JF, the stress-energy tensor is covariantly conserved (Faraoni 2004). We assume the standard equation of state 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 with ρ 0 the mean density of the component at the present time, t 0 , when a(t 0 ) = a 0 = 1.

A spatially flat universe: Analytic solution
Here, we solve the field equations, Eqs. (32)-(34), in the special case of a spatially flat universe with k = 0. We assume that the universe only contains baryonic matter with density ρ b and negligible pressure, p = 0, namely w = 0 in the equation of state, Eq. (35). After substituting Eq. (32) into Eq. (33) and dividing by The solution of the above system of coupled equations together with the equation of the mass-energy conservation, Eq. (36), determines the time evolution of both the scalar field and the scale factor. We simplify the field equations by introducing the modified cosmological parameter where Ω b = ρ b /ρ cr ≡ ρ b (8πG/3H 2 ) is the density parameter associated with the homogeneous baryonic density ρ b . The density parameter Ω is analogous to the total matter density parameter of 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/ϕ, unlike the standard model, where the baryonic and the dark matter components simply add up. Similarly, we define The parameter Ω Ξ is analogous to the standard cosmological parameter Ω Λ = Λ/3H 2 , with Λ the cosmological constant in the GR field equations; therefore, in CRG, Ξ exactly plays the role of Λ in the standard model. At t = t 0 , the values of the two parameters Ω and Ω Ξ are Ω 0 = 16πGρ b0 /3H 2 0 ϕ 0 and Ω Ξ0 = Ξ/3H 2 0 . We further introduce two deceleration parameters, q and q ϕ , related to the scale factor and the scalar field, respectively: With the above definitions, the field equations, Eqs. (37)-(39), become It is convenient to define the new quantity ζ ≡φ Hϕ .
The field equations, Eqs. (43)-(45), now become Combining the first two equations yields namely, with the definitions of q and ζ, This equation can be integrated to obtain ϕ = C/(a 2ȧ ), with C an integration constant. Its value can be found from the boundary conditions a(t 0 ) = 1, H(t 0 ) = H 0 , and ϕ(t 0 ) = ϕ 0 . We obtain C = H 0 ϕ 0 and thus Moreover, since 3Ha 3 = da 3 /dt, we can write the scalar field as This expression shows that the scalar field ϕ is inversely proportional to the rate of the variation of the volume of the universe. As shown in Appendix C, adding up Eqs. (47) and (49) yields two solutions for dH/dt: each solution being a first order differential equation for H. Hereafter, we refer to the solution corresponding to the upper (+) sign as CRG+, and to the lower (−) sign as CRG-, respectively. The integration of Eq. (54) yields with the upper and lower signs corresponding to CRG+ and CRG-, respectively. The integration constant, C 1 , can be found from the boundary condition at t = t 0 : Furthermore, the scale factor can be derived by integrating Eq. (55) and C 2 is an integration constant, whose value can be fixed by the condition a 0 = 1: with Inverting Eq. (60) yields which implies that the cosmological parameters Ω 0 and Ω Ξ0 satisfy the relation The parameter Ω 0 must be positive, because it is proportional to the mass-energy density ρ b0 , according to Eq. (40); therefore, the only physically viable solutions of these inequalities are The parameter Ξ must thus be positive. The first test of the viability of CRG is its ability to properly describe the Hubble diagram of the observed SNeIa at high redshift. Using the scale factor in Eq. (57), the luminosity distance can be computed exactly, for both signs appearing in Eq. (57) (see Appendix D). As in the standard cosmological model, the values of the cosmological parameters can be set by modelling the SNIa Hubble diagram (see Appendix E).

The equation of state of the effective dark energy
Equation (13) is the general expression of the stress-energy tensor of the scalar field. By using the FLRW metric of Eq. (31), Eq. (13) becomes If we consider T ϕ µν to be in the form of a perfect fluid, T ϕ µν = ρ ϕ + p ϕ u µ u ν + p ϕ g µν , the 00-and i j-components of Eq. (64) lead to the effective density and pressure associated with the scalar field: In order to derive the density and pressure associated with the effective dark energy in CRG, we can combine Eq. (32) with Eq. (34), and Eq. (33) with Eq. (34), to rewrite the modified Friedmann Eqs. (32)-(33) as By comparing these equations with those of a general TeS theory with a non-minimal coupling between the scalar field and the metric, where ϕ 0 is the value of the scalar field computed at the present epoch, 11 we derive the general expressions where ρ = ρ m +ρ r and p = p r = ρ r /3 are the density and pressure associated with matter, ρ m , and radiation, ρ r and p r , respectively (Frusciante & Perenon 2020).
Based on these definitions, we can evaluate the parameter associated with the equation of state of the effective dark energy, w DE = p DE /ρ DE . Its analytical expression is reported in Eq. (F.4) of Appendix F. At the present epoch, for H = H 0 and a = a 0 = 1, we obtain where the upper and lower signs refer to the CRG+ and the CRG-solutions, respectively. The value which best fits the observational data is w DE = −1 (Aghanim et al. 2020c), which is consistent with the cosmological constant of the ΛCDM model. By assuming Ω 0 ∼ 0.3, we find Ω Ξ0 ∼ 0.64, which is only admitted by CRG-and is near the upper bound of Eq. (63).
Values of w DE slightly different from −1 can nonetheless be accommodated by the data (e.g. Amendola & Tsujikawa 2015;Copeland et al. 2006;Frusciante & Perenon 2020;Wen et al. 2018;Capozziello et al. 2006;Gerardi et al. 2019). Observational constraints on the equation of state of the effective dark energy generally depend on the model used to describe its effects (but see Gerardi et al. 2019, for a model-independent reconstruction). We can gain some insight by adopting the parametrization w DE (z) = w 0 + w a z/(1 + z) (Chevallier & Polarski 2001;Linder 2003). At the present epoch z = 0, this parametrization only depends on w 0 , and it is therefore sufficiently general to account for a broad range of dark energy models, either with w DE < −1 (phantom models) or w These values of Ω Ξ0 close to the upper limit (1 + Ω 0 )/2 (Eq. 63) are consistent with the constraints we obtain from the SNIa data that we investigate in Appendix E. Future observations aimed at investigating the large-scale structure of the Universe and the evolution of dark energy (Amendola et al. 2018) will further constrain w 0 , and hence the value of Ξ.

Additional remarks on the evolution of the scalar field
By using Eqs. (40), (41), and (36) for a matter-dominated universe (p = 0), we can write the field equation of the scalar field, Eq. (52), as The right-hand side of this equation is negative at the present time t 0 , for Ω 0 = 0.3 and Ω Ξ0 = 0.65. In other words, the contribution of Ω Ξ0 dominates over Ω 0 . However, the Ξ term plays a decreasingly important role at increasing redshifts with the righthand side of the equation above becoming positive, by adopting H(t) from Eq. (55), at times t earlier than In the radiation-dominated era, however, the right-hand side of Eq. (34) is zero, and therefore Ξ may become an important source term in the scalar field equation.
Big Bang nucleosynthesis could potentially constrain the value of Ξ and, more generally, the scalar potential. Scalar-tensor theories, and specifically Brans-Dicke-like models, provide a time-varying effective gravitational constant. This feature alters the rate of expansion of the universe and therefore its temperature, which both contribute to regulate the primordial nucleosynthesis. In particular, the observed abundances of 2 H, 3 He, and 7 Li set a bound on the present matter density, whereas the abundance of 4 He constrains the rate of change of ϕ (Arai et al. 1987). For TeS theories without a self-interaction potential, there is an attractor mechanism towards GR (Damour & Pichon 1999) as the scalar field, and thus G, remains approximately constant during the radiation-dominated epoch. Therefore, in these models, the nucleosynthesis of light elements occurs similarly to the standard cosmological scenario, but with a different expansion rate (Damour & Gundlach 1991;Casas et al. 1992a,b;Serna et al. 1992;Clifton et al. 2005). The primordial nucleosynthesis in the presence of a self-interaction potential -as in CRG -has been extensively studied in the literature (see, e.g., Arai et al. 1987;Uzan 2003;Larena et al. 2007;Coc et al. 2006;Iocco et al. 2009;Clifton et al. 2012, and references therein), also as a possible solution to the problem of the 7 Li abundance. An extensive analysis of whether CRG is consistent with the observed abundances of light elements is left to future work.

Conclusion
Refracted gravity 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 (ρ), a monotonic function of the mass density ρ (Matsakos & Diaferio 2016). Here, we propose a covariant extension of RG, CRG, and show that it 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 V(ϕ) of the scalar field ϕ and by the general differentiable function W(ϕ) appearing in the Lagrangian density. For CRG, we adopt V(ϕ) = −Ξϕ, with Ξ a normalization constant, and W(ϕ) = −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 ϕ = 2 .
When used to describe an expanding and homogeneous universe, the scalar field ϕ is also responsible for the observed accelerated expansion of the Universe. The cosmological density associated with the cosmological constant in the standard model, Ω Λ0 = Λ/3H 2 0 , is now replaced by the CRG parameter Ω Ξ0 = Ξ/3H 2 0 . It follows that the normalization constant Ξ plays the role of the cosmological constant Λ of the standard model. Moreover, the Hubble diagram of high-redshift SNeIa and current observational constraints on the parameter w DE of the equation of state of the effective dark energy, that we derive from the stress-energy tensor associated with the scalar field, suggest that, in a universe with a flat geometry, Ω Ξ0 ∼ (1 + Ω 0 )/2, with Ω 0 the ratio between the baryonic matter and the homogeneous scalar field at the present time. The parameter Ω 0 can be identified with Ω 0 = 0.31 of the standard model (Aghanim et al. 2020c). We find Ω Ξ0 = 0.650 +0.005 −0.085 at the 90% confidence level for the CRG+ solution. It thus follows that Ξ and Λ have comparable values.
In addition, in the weak-field limit, Ξ sets an acceleration scale, (2Ξ) 1/2 ∼ 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 et al. 2022). This acceleration scale is indeed present in real systems (Chae et al. 2020;McGaugh et al. 2018) and is comparable to the acceleration a 0 introduced in MOND (Milgrom 1983a,b,c). Therefore, being Ξ ∼ Λ, the known relation a 0 ∼ Λ 1/2 (Milgrom 1999(Milgrom , 2020 naturally appears in CRG. CRG provides a connection between phenomena generally attributed to dark matter and dark energy separately and it thus belongs to the family of modified gravity models that connect the two dark sectors within a unified scenario. The same property is currently shared by other models, including some f (R) theories (Sotiriou & Faraoni 2010), "quartessence" theories  or generalized Chaplygin gas models (Bento et al. 2002;Zhang et al. 2006).
Assessing if the evolution of the universe described by CRG is consistent with observations requires extensive additional work. For example, the power spectrum of the temperature anisotropies of the CMB and the power spectrum of matter density perturbations are required to constrain the Lagrangian density and its parameters (see, e.g., Noller & Nicola 2019b;Huterer et al. 2015). Similarly, the analysis of the evolution of the density perturbations, at least in the linear perturbation theory (Di Porto & Amendola 2008;Bueno Sanchez & Perivolaropoulos 2011;Pace et al. 2014;Kofinas & Lima 2017), and of the role of the scalar field and its perturbations on structure formation are of crucial importance to test the viability of CRG.
Finally, since the scalar field drives the accelerated expansion of the Universe, we need to further investigate the connection of CRG with current dark energy models (Pettorino et al. 2005;Capozziello et al. 2006;Frusciante & Perenon 2020). Specifically, we find that the parameter w DE of the equation of state of the effective dark energy depends on redshift, unlike the standard cosmological model. When tested with measures from the upcoming Euclid mission (Amendola et al. 2018), these predictions could discriminate CRG from current dark energy models. We plan to tackle all these issues in future work.
The ∇ µ ∇ ν ϕ components are For a static non-relativistic fluid, whose pressure p is negligible compared to its density ρ, p ρ, the equation u α u α = −1 implies u 0 = (−g 00 ) 1/2 ; therefore, the components of the energymomentum tensor T µν = ρu µ u ν are and its trace is Inserting these relations into Eqs. (10) and (11), we obtain the system of equations ∇ · (ϕ∇Φ) 8πGρ , (A.16) where ρ(r) = ρ s (r) + ρ bg , with ρ s (r) being the density profile of the source in a homogeneous background with density ρ bg , and m(< r) is the enclosed mass within radius r. Newtonian gravity is recovered when ϕ = 2 (i.e. = 1). For spherically-symmetric systems, Eqs. (A.17) and (A.18) are less straightforward to solve than Eq. (A.16). In the following, we make a few approximations to derive the generic behaviour of the gravitational field dΦ/dr and the scalar field ϕ in the two limits of small and large distances from the spherical source.
A.1.1. The gravitational field at small distances from the spherical source In the Newtonian regime where ρ s ρ bg , we may approximate ϕ 2−ϕ 1 , where |ϕ 1 | 1. By ignoring all the terms of the order O(ϕ 2 1 ), O(ϕ 1 )O(Φ), and O(ϕ 1 )O(U), Eq. (A.18) yields where we also ignored the terms containing Ξ: as we show in Sect. 4, we can identify Ξ 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. In the vacuum ρ bg = 0, and the RG field reduces to the Newtonian field for ϕ = 2, in agreement with Eq. (A.16), and it is larger than the Newtonian field for ϕ < 2.
A.1.2. The gravitational field at large distances from the 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 ρ bg approximates the mean density of the Universe, assuming that the source is isolated. We can thus treat the contribution of the source, ρ s , as a small perturbation to the total density field ρ = ρ bg + ρ s , with ρ s ρ bg . We can similarly assume that the scalar field reaches a mean cosmic value ϕ bg = const and write the scalar field as ϕ = ϕ bg + ϕ 1 , with ϕ 1 ϕ bg . In the limit ρ s → 0 and ϕ 1 → 0, Eqs. (A.17) and (A.18) simplify to We can derive how the scalar field approaches the cosmic constant value ϕ bg by assuming that Eq. (A.24) already holds at distances where ϕ is not yet asymptotically constant, i.e. ∇ 2 U 2πGρ/ϕ. By adding up Eqs. (A.17) and (A.18), and using Eq. (A.16), we obtain In spherical symmetry, we can recast this equation as 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ϕ/dr < 0. We now explore the two limits of the large and small gravitational field dΦ/dr at large distances from the source. In the limit dΦ/dr (2Ξ − 8πGρ/ϕ) 1/2 , Eq. (A.27) reduces to dϕ/dr ∝ −2ϕ dΦ/dr, a dependence broadly approaching the result derived in the previous section for small distances from the source (Eqs. A.19 and A.21).