Refracted Gravity Solutions from Small to Large Scales

If visible matter alone is present in the Universe, general relativity (GR) and its Newtonian weak field limit (WFL) cannot explain several pieces of evidence, from the largest to the smallest scales. The most investigated solution is the cosmological model $\Lambda$ cold dark matter ($\Lambda$CDM), where GR is valid and two dark components are introduced, dark energy (DE) and dark matter (DM), to explain the $\sim$70\% and $\sim$25\% of the mass-energy budget of the Universe, respectively. An alternative approach is provided by modified gravity theories, where a departure of the gravity law from $\Lambda$CDM is assumed, and no dark components are included. This work presents refracted gravity (RG), a modified theory of gravity formulated in a classical way where the presence of DM is mimicked by a gravitational permittivity $\epsilon(\rho)$ monotonically increasing with the local mass density $\rho$, which causes the field lines to be refracted in small density environments. Specifically, the flatter the system the stronger the refraction effect and thus, the larger the mass discrepancy if interpreted in Newtonian gravity. RG presented several encouraging results in modelling the dynamics of disk and elliptical galaxies and the temperature profiles of the hot X-ray emitting gas in galaxy clusters and a covariant extension of the theory seems to be promising.

The most investigated solution to explain this phenomenology is to assume the presence of a nonbaryonic and cold, i.e., nonrelativistic at the epoch of decoupling from radiation, form of matter which only gravitationally interacts with baryonic matter: the cold dark matter (CDM) [12].However, the only presence of dark + baryonic matter is again not sufficient to account for all the pieces of evidence in the Universe, which suggest a further discrepancy of ∼70%.The most studied cosmological model is the ΛCDM, which assumes GR as the gravity theory and introduces two dark constituents besides visible matter: dark energy (DE) and dark matter (DM), which can explain the ∼70% and ∼25% of the mass-energy budget of the Universe.Specifically, DE is an exotic fluid with negative pressure which justifies the accelerated expansion of the Universe, as observed from the Hubble diagram of Ia Supernovae (SNe Ia) [13,14], and it can be identified with the cosmological constant Λ.
Even if ΛCDM can account for the majority of the observations in the Universe, it presents some problems, both on large and small scales.On large scales, we observe the cosmological constant problem [15,16], the coincidence problem [17], and the tensions between the values of some cosmological parameters measured from probes of the late and early Universe [1,18,19].On small scales, we observe several discrepancies between ΛCDM simulations and observations.These generate the cusp/core, missing satellites, too-big-to-fail, and planes of satellite galaxies' problems (e.g., see [20,21] for a review).On a galaxy scale, another issue related to the dynamics of galaxies is present in ΛCDM.This is the "disk-halo conspiracy", since the relative contributions of the baryonic stellar disk and the DM halo to the overall rotation curves of disk galaxies are degenerate to each other, that is, their contribution is not univocal.Several pieces of evidence suggest the "maximum-disk hypothesis" to solve this issue [22][23][24][25][26][27], which consists in assuming that the baryonic stellar disk maximally contributes to the innermost part of the rotation curve.Yet, this assumption is only valid for high-surface-brightness (HSB) galaxies, whose rotation curves steeply rise in their inner regions and are modelled with a cuspy DM density profile [28,29].On the contrary, low-surface-brightness (LSB) galaxies, generally dwarf and dwarf spheroidal (dSph) galaxies, are ones of the darkest objects at a galaxy scale in the Universe and they are DM-dominated as a result, even in their central regions: they have an inner velocity dispersion σ ∼ 10 km s −1 , an order of magnitude larger than the velocity dispersion σ ∼ 1 km s −1 for systems with the same luminosity and scale radius (∼100 pc) at equilibrium [30].Their rotation curves, slowly rising toward their flat region, are described by a cored DM density profile and not by a maximum-disk model (e.g., [31]), which does not naturally emerge in ΛCDM cosmological simulations and gives rise to the cusp/core problem (e.g., [20,21]).Instead, there are systems with baryonic masses similar to those of dwarf galaxies, the globular clusters (GCs), which are nearly DM-free [32][33][34][35][36][37][38].
It is remarkable to cite also the presence of some observed regularities on a galaxy scale, which are hard to explain in a ΛCDM paradigm, where a stochastic merging process of structures' formation is invoked.Some of these regularities are very tight scaling relations between a property of the DM and of the baryonic matter in galaxies, which might be counter-intuitive since DM represents ∼90% of the galaxy content, whereas baryonic matter only represents ∼10%.Among these relations, we can mention (1) the baryonic Tully-Fisher relation (BTFR) [39], (2) the mass discrepancy-acceleration relation (MDAR) [40], and (3) the radial acceleration relation (RAR) [41].The three relations see the emergence of the same acceleration scale a 0 ≃ 1.2 × 10 −10 m s −2 , which, if expressed in natural units, is a 0 ∼ H 0 ∼ √ Λ, suggesting a connection between the DM and DE sectors.These coincidences are even less intuitive in a ΛCDM context.
An alternative way to explain these discrepancies or regularities is provided by alternative theories of gravity, without the addition of any dark constituent.One of the most investigated modified theories of gravity formulated in a nonrelativistic way is modified Newtonian dynamics (MOND) [42][43][44].MOND assumes a modification of the law of gravity dependent on the value of the background acceleration: when its value goes below the acceleration scale a 0 , the gravitational field departs from Newtonian gravity, being subjected to a boost that mimics the effect of DM.MOND not only explains more intuitively but actually predicts several aspects of the dynamics of galaxies, such as the flatness of the rotation curves of disk galaxies and the three mentioned scaling relations.Matsakos and Diaferio [45] proposed in 2016 a different approach.They formulated refracted gravity (RG), a classical modified theory of gravity which does not assume the presence of DM and is regulated by the value of the local mass density ρ, rather than of the acceleration a.The Poisson equation of RG is modified at first member by the presence of the gravitational permittivity, ϵ(ρ), a monotonic increasing function of ρ which boosts the gravitational field in regions where the density goes below a critical value, reproducing the effect of DM in Newtonian gravity.A covariant version of RG was recently formulated [46] and it seems to describe both the DE and DM sectors with a single scalar field and to reproduce the Hubble diagram of SNe Ia.
RG has obtained some encouraging results in modelling the dynamics of galaxies and galaxy clusters, as well as on the covariant side.The paper develops as follows.Section 2 describes RG theory.Sections 3-6 recap the main analyses and results obtained with RG in the field of disk galaxies, elliptical galaxies, galaxy clusters, and covariant RG, respectively.Section 7 discusses the future projects of RG and concludes the paper.

Refracted Gravity
RG is a theory of modified gravity formulated in a classical way, which can be interpreted in analogy to electrodynamics in matter [45]: as an electric field line suffers a change both in direction and in magnitude when it crosses a dielectric medium with a nonuniform permittivity, a gravitational field line suffers the same changes when it passes from a high-density to a low-density environment.This behaviour of the gravitational field is encoded in this modified Poisson equation (Equation (2.3) in [45]): where ϕ is the gravitational potential, and ϵ(ρ) is the gravitational permittivity.The following asymptotic limits are adopted for the gravitational permittivity (Equation (2.5) in [45]): where 0 < ϵ 0 ≤ 1, and ρ c are the gravitational permittivity in vacuum and the critical density, respectively, two of the three free parameters of the theory.When ρ ≫ ρ c , the Newtonian Poisson equation is recovered (Equation (2.7) in [45]): whereas, when ρ ≪ ρ c , we are in a fully RG regime.
The RG gravitational field has a different behaviour for spherical and flat systems.For spherical systems, the RG Poisson equation (Equation ( 1)) reduces to (Equation (3) in [47]): where M(< r) is the mass of the spherical system enclosed within the spherical radius r, and ∂ϕ N /∂r is Newtonian gravitational field.In this case, the RG field direction and r-dependence remain Newtonian, and the RG field magnitude increases compared to the Newtonian one where ρ ≪ ρ c , i.e., where ϵ(ρ) → ϵ 0 (right panels of Figure 1).If we consider a solid sphere immersed in vacuum (ρ = 0) with radius r s and uniform density ρ s ≫ ρ c , the RG gravitational field inside and outside the sphere is given by [45]: where g = ∂ϕ/∂R and g N = ∂ϕ N /∂R are the RG and Newtonian gravitational fields, respectively.The gravitational field has the Newtonian direction everywhere but is enhanced by 1/ϵ 0 over the Newtonian field outside the sphere.The analogy with electrodynamics in matter, i.e., the refraction of the field lines where the density decreases, is observed for flattened systems.Expanding the left-hand side of Equation (1), we obtain (Equation (4) in [48]): In this configuration, the RG field depends both on the density ρ (second term on the left-hand side of Equation ( 6)) as for spherical systems, and on its gradient (first term on the left-hand side of Equation ( 6)).In particular, the term " ∂ϵ ∂ρ ∇ρ • ∇ϕ" is different from zero in nonspherical configurations and causes the field lines to refract.The acceleration boost in the external regions of disk galaxies that in Newtonian gravity is explained by the presence of DM, is explained in RG by the refraction of the field lines toward the equatorial plane of the disk caused by the low-density regions above and below the disk plane (left panels of Figure 1).Following Equation ( 6), RG predicts that the flatter the system, the larger the mass discrepancy if interpreted in Newtonian gravity.This might explain a positive correlation between ellipticity and DM content of elliptical galaxies [49,50] and the different DM quantity of GCs and dwarf galaxies, where the former are nearly spherical and DM-free (e.g., [35]) and the latter are flatter and one of the darkest objects in the Universe (e.g., [30]) (see Section 1).It has to be specified that the formulae reported for flattened systems in the left panels of Figure 1 are equivalent formulae for spherical systems to approximately illustrate the dependence of the field on the distance from the galaxy centre.Where ρ ≪ ρ c in flat systems, e.g., in the outskirts of disk galaxies, the radial gravitational field assumes the asymptotic trend (Equation (B.23) in [48]): where a 0 is the MOND critical acceleration set by the normalisation of the BTFR [51]: where v f is the asymptotic value of the flat part of the rotation curve of disk galaxies and M is the total baryonic mass of the disk galaxy.Specifically, Matsakos and Diaferio [45] derived the relation v 4 f = GbM (Matsakos and Diaferio [45] used the lowercase font for the baryonic mass m, whereas I adopt the uppercase font for the same symbol across the paper), where they put the acceleration parameter b in relation with the MOND acceleration scale a 0 to be consistent with the normalisation of the observed BTFR [51,52].In this regime, the RG field trend deviates from the Newtonian inverse square law.This asymptotic limit coincides with the MOND asymptotic limit for the gravitational field where a ≪ a 0 , which might indicate that RG shares the majority of MOND's successes on a galaxy scale.The difference between RG and MOND is better observed in spherical systems, where the field trend in RG is Newtonian (∝ R −2 ) even where ρ ≪ ρ c , preserving the Gauss theorem, but the field trend in MOND is ∝ R −1 where a ≪ a 0 , as in flattened systems.
In the regions of flat systems where ρ ≪ ρ c , besides the asymptotic limit given by Equation (7) for the radial gravitational field, the following asymptotic limit for the vertical gravitational field holds (Equation (B.25) in [48]): which means that the field is redirected parallel to the equatorial plane of the flat system at the height where condition ( 9) is reached.If we consider a solid disk immersed in vacuum (ρ = 0) with radius R d , height h d from its equatorial plane (2h d of total height), and uniform density ρ d ≫ ρ c , its RG gravitational field remains Newtonian within the disk and reaches the asymptotic limits ( 7) and ( 9) outside the disk.Specifically, condition ( 9) is accomplished at distance z = ±h d from the equatorial plane of the disk.Different functions can be adopted for the gravitational permittivity, such that it accomplishes the asymptotic limits of Equation (2).For all the analyses performed with RG, the following smooth step function of the mass density ρ was employed (Equation (4.1) in [45]), Equation ( 5) in [48], and Equation ( 4) in [47]: represented in Figure 2.Besides ϵ 0 and ρ c , a third free parameter, Q, is present in RG, where Q regulates the steepness of the transition between the two asymptotic limits of Equation ( 2).In Figure 2, Equation ( 10) is represented for three values of Q, showing that the larger its value, the steeper the transition between the two asymptotic regimes of Equation (2).The three RG free parameters are supposed to be universal.

Dynamics of Disk Galaxies
Matsakos and Diaferio [45] presented some preliminary encouraging results, according to which RG properly models the BTFR (Equation ( 8)) and the MDAR of galaxies, and the rotation curves of the HSB disk galaxy NGC 6946 and LSB disk galaxy NGC 1560.Cesare et al. [48] demonstrated that RG can properly model the dynamics of disk galaxies.For their analysis, they considered the rotation curves and the vertical velocity dispersion's radial profiles of 30 disk galaxies in the DiskMass Survey (DMS) [53] at low redshift.They chose this sample since these galaxies are close to being face-on, which allowed the measurement of both their rotation curves and vertical velocity dispersions.In this way, Cesare et al. [48] could obtain a single constraint of the RG parameters from two kinematic profiles taken at the same time rather than from the rotation curves alone, which provided a more robust analysis.These results are also summarised in [54].
Cesare et al. [48] modelled with RG (1) the rotation curves alone of each DMS galaxy, (2) the rotation curves and vertical velocity dispersions at the same time of each DMS galaxy, (3) the rotation curves and vertical velocity dispersions of all 30 DMS galaxies at the same time, and (4) the RAR of DMS galaxies.

Mass Model
Cesare et al. [48] modelled the baryonic mass density profile ρ, which generates the total gravitational potential ϕ in the hypothesis that DM is not present, with (1) an axisymmetric exponential disk for the stars, (2) a spherical stellar bulge, and (3) two razorthin disks for the atomic and molecular gas.The 3D mass density of the stellar disk was modelled with an exponential profile of this kind (Equation (7) in [54]): where I d (R) is the surface brightness radial profile of the stellar disk, R and z are the cylindrical coordinates (R is the radius projected on the sky, oriented along the major axis of the disk, and z is the vertical coordinate, oriented perpendicular to the disk equatorial plane), and Υ and h z are the stellar mass-to-light ratio and the disk scale height, two of the five free parameters of the dynamical model adopted by [48].
The surface brightness of the disk was modelled with a linear interpolation of the measured surface brightness data to catch the specific features of the luminous matter distribution, which generally correspond to features in the rotation curve, following "Renzo's rule" [55].The disk surface brightness in the innermost galaxy regions, where the bulge is dominant, was modelled with an exponential profile [56,57] (Equation (A.1) in [48]): obtained by fitting the surface brightness points of the outermost stellar disk, where the stellar disk alone dominates.DMS galaxies are disk-dominated, the average stellar bulge-over-total luminosity ratio in the K-band being around 0.09.For this reason, the bulge was approximated as a sphere, introducing negligible systematic errors.The 3D mass density of the bulge was modelled with this Abel integral (Equation (9) in [54]): where I b (R) is the surface brightness radial profile of the bulge, and is the 3D radius.The bulge contribution being subdominant, Cesare et al. [48] assumed the mass-to-light ratio of the bulge equal to the mass-to-light ratio of the disk, Υ, without introducing an additional degree of freedom.They modelled the surface brightness of the bulge with a Sérsic spherical profile [58] (Equation (A.2) in [48]): To account for the seeing, which is not negligible since the measurements were taken from the ground with the 3.5 m diameter Calar Alto Observatory, they modelled the observed surface brightness of the bulge with Equation ( 15) convolved with a Gaussian point spread function.
The 3D mass density of the atomic and molecular gas were modelled with two razorthin disks (Equation (A.14) in [48]): and where Σ atom (R) and Σ mol (R) are the surface mass densities of the atomic and molecular gas, respectively, which Cesare et al. [48] modelled with a linear interpolation of the data as for the stellar disk, to catch the features of the rotation curve, and δ(z) is the Dirac δ function.

Dynamical Model
From the total 3D mass density profile given by Equation (18), Cesare et al. [48] derived the RG gravitational potential by solving the RG Poisson Equation (1).To solve the RG Poisson equation, they used an iterative Poisson solver based on the successive-overrelaxation method, given by combining the Jacobi and Gauss-Seidel methods [59].From the RG potential, they derived the models for the rotation curve, v(R) (Equation ( 6) in [48]), and the vertical velocity dispersion radial profile, σ z (R) (Equation ( 13) in [48]), where Equation ( 20) is the vertical velocity dispersion for a vertically decreasing exponential disk whose density follows Equation (11).
The dynamical model has five free parameters: the stellar mass-to-light ratio Υ, the stellar disk scale height h z , and the three RG parameters, ϵ 0 , Q, and ρ c .They estimated the free parameters of the model firstly from the rotation curve alone of each galaxy and then from the rotation curve and the vertical velocity dispersion taken at the same time of each galaxy with a Monte Carlo Markov Chain (MCMC) algorithm with a Metropolis-Hastings acceptance criterion.They adopted the following priors for the free parameters of the model.For Υ, they adopted a Gaussian prior based on the stellar population synthesis (SPS) model of Bell and de Jong [60], keeping into account that the surface brightness was measured in the K-band, where the Gaussian tail for Υ < 0 was set to zero.For h z , they adopted a Gaussian prior peaking at h z,SR , a scale height obtained with the following relation, derived from a combined sample of 60 edge-on late-type galaxies by Bershady, et al. [61] (Equation (1) in [61] and Equation (A.12) in [48]): where the term ±0.095 indicates the ∼25% intrinsic scatter.Estimating h z,SR was essential to have a comparison reference for DMS galaxies where, due to their nearly face-on configuration, their h z could not be directly measured.For the three RG parameters, they adopted a flat prior in the following intervals: [0.10, 1.00] for ϵ 0 , [0.01, 2.00] for Q, and [−27.00,−23.00] for log 10 [ρ c (g/cm 3 )].

Results
RG modelled the kinematic profiles of DMS galaxies with sensible parameters.Figure 3 shows the RG models (blue solid lines) of the rotation curves (panels (a) and (b)) and the vertical velocity dispersion profiles (panels (c)) against the DMS data (red dots with error bars) for three DMS galaxies (UGC 3091, UGC 3701, and UGC 4256).The blue curves for the rotation curves in panels (a) were computed with the parameters estimated from the rotation curves alone of each galaxy.The blue curves for the rotation curves in panels (b) were computed with the parameters found by modelling the rotation curves and vertical velocity dispersion profiles of each galaxy at the same time.With the same parameters, the vertical velocity dispersion profiles were computed in panels (c).In panels (c), the cyan solid lines are the RG vertical velocity dispersions calculated with the same parameters as the blue curves but with disk-scale heights h z set to h z,SR (Equation ( 21)).In all panels, the dashed magenta vertical lines identify the bulge effective radius (R e in Equation ( 15)), and the dashed green vertical lines represent the bulge radius adopted in the disk-bulge decomposition for the surface brightness modelling (see Section 3.1).We can see that RG properly modelled both kinematic profiles and that the substructures of the rotation curves were properly reproduced following Renzo's rule (see Section 3.1).
Both in the only-rotation curve and in the combined rotation curve + vertical velocity dispersion analyses, Cesare et al. [48] found mass-to-light ratios Υ generally consistent with the SPS model of Bell and de Jong [60].In the only-rotation curve analysis, the estimated h z 's were generally consistent with the h z,SR derived from Equation (21).Instead, in the combined analysis, the estimated h z 's were systematically smaller than the correspondent h z,SR .Angus et al. [62] found a similar result by modelling the DMS galaxies with QUMOND [63], a modified gravity version of MOND theory.Yet, this result seems to not indicate a flaw in the two modified theories of gravity but instead, an observational bias.Indeed, the disk scale height h z found by Cesare et al. [48] and by Angus et al. [62] was estimated from the vertical velocity dispersions, which are measured with spectroscopy from a signal derived from a mixed stellar population of young and old stars, where the signal from the young population is dominant [64].The young stellar population distributes in a disk thinner and with a smaller vertical velocity dispersion σ z than the old stellar population [64].Instead, the disk scale heights of edge-on galaxies, from which Equation ( 21) is derived to determine h z,SR , are directly measured from near-infrared photometry, whose signal is dominated by the old stellar population [64], having higher h z and σ z than the young stellar population.It is possible to simulate the vertical velocity dispersion of the old stellar population by artificially increasing the vertical velocity dispersion of the young stellar population by a proper factor.Cesare et al. [48] computed this factor, which was equal to 1.63, in agreement with the work of Aniyan et al. [64].Cesare et al. [48] artificially increased this factor to determine the vertical velocity dispersions of five DMS galaxies, obtaining disk scale heights in agreement with the correspondent h z,SR .
UGC 4380   The cyan solid lines are the RG vertical velocity dispersion profiles calculated with the same parameters as the blue curves but with disk-scale heights h z = h z,SR (Equation ( 21)), and the dashed magenta and green vertical lines show the bulge effective radius (R e in Equation ( 15)) and the bulge radius adopted in the disk-bulge decomposition for the surface brightness modelling, respectively.The figure is readapted from Figures D.1 and D.2 in [48].Credit: Cesare V., Diaferio A., Matsakos T., and Angus G., A&A, 637, A70, 2020, reproduced with permission © ESO.
The three RG parameters derived from the rotation curves and the vertical velocity dispersions of each galaxy are in agreement with each other: their average errors are larger than the standard deviation of their distributions.This suggests their universality, as should be from the theory formulation.To better test their universality, Cesare et al. [48] estimated a unique combination of RG free parameters from the two kinematic profiles of all 30 galaxies taken at the same time, to verify their agreement with the average of the RG parameters found from the rotation curves and the vertical velocity dispersions of each galaxy.The mean RG parameters derived from the single galaxies with their mean errors were ϵ 0,Mean,DMS = 0.56 ± 0.16, Q Mean,DMS = 0.92 ± 0.71, log 10 [ρ c (g/cm 3 )] Mean,DMS = −25.30± 1.22 (purple squares with error bars in Figure 4), whereas the unique combination of free parameters was ϵ 0,Unique,DMS = 0.661 +0.007 −0.007 , Q Unique,DMS = 1.79 +0.14 −0.26 , log 10 [ρ c (g/cm 3 )] Unique,DMS = −24.54+0.08 −0.07 , where their estimates and error bars were the medians (green dots in Figure 4) and the 15.9 and 84.1 percentiles of their posterior distributions (blue shaded areas in Figure 4).In Figure 4, the yellow, red, and black contours show the 1σ, 2σ, and 3σ levels.The code used to find this unique combination of RG parameters was written in C++ and parallelised with OpenMP.It is publicly available on GitHub under the name "astroMP" (https://github.com/alpha-unito/astroMP,accessed on 7 November 2019) and it is described in [65,66].The unique combination of RG parameters was found with an approximate procedure, in which the mass-to-light ratios and the disk scale heights were kept fixed to the values found from the single galaxies.The fact that this unique combination of free parameters provided a limited worsening of the agreement between the RG models and the kinematic data and that the unique ϵ 0 , Q, and log 10 [ρ c (g/cm 3 )] agreed with the correspondent RG parameters estimated from the single galaxies within 0.63σ, 1.19σ, and 0.62σ, respectively, further enforced the hypothesis of the universality of these parameters.With a proper exploration of the 63-dimensional parameter space for this galaxy sample (one Υ and one h z for each of the 30 galaxies and three RG parameters for all the galaxies), a better consistency between model and measurements and a more precise estimate of the unique combination of RG parameters might be achieved.
With the parameters obtained from the two kinematic profiles of each galaxy, Cesare et al. [48] also modelled the RAR of the DMS galaxies.McGaugh et al. [41] fitted the RAR of 153 edge-on disk galaxies in the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample [67] with this relation (Equation ( 4) in [41] and Equation (20) in [48]): where g obs (R) is the observed centripetal acceleration derived from the measured rotation curve v obs (R) (Equation ( 19) in [48]), g bar (R) is the Newtonian acceleration due to baryonic matter alone, and g † is the only free parameter of the model.This relation properly interpolated the SPARC data with g † = [1.20 ± 0.02 (random) ± 0.24 (systematic)] × 10 −10 m s −2 , which agrees with the MOND acceleration scale a 0 within 1σ.
The RG models of the RAR for the DMS galaxies (blue solid lines in Figure 5) properly interpolated the measured RAR of the DMS data (red dots with error bars in Figure 5) and reproduced the asymptotic limits of Equation (22).However, RG reproduced the RAR of DMS galaxies with a too large intrinsic scatter (0.11 dex vs. 0.057 dex as found by Li et al. [68] from SPARC galaxies) and with residuals strongly correlated, largely at more than 5σ, with some galaxy properties dependent on the distance from the galaxies' centre (the radius, the stellar surface density profile, and the gas fraction profile), apparently at odds with the results found by Lelli et al. [69] for SPARC galaxies.Yet, the RAR of DMS galaxies also presented some, although weaker, correlations between their residuals from Equation (22) and some galaxy properties.A galaxy sample such as the DMS, made of only 30 galaxies which are also close to being face-on and consequently, with not very accurately measured rotation curves, might not be particularly suitable to investigate the RAR, but a larger sample of edge-on disk galaxies, such as SPARC, should be considered for further analyses.However, further investigation has to be conducted in this sense, since this issue might not depend on the galaxy sample but on the RG theory itself.

Dynamics of Elliptical Galaxies
In addition to the dynamics of flat systems (disk galaxies), RG has also been tested on the dynamics of spherical systems.To test RG for this class of objects, Cesare et al. [47] considered a sample of three nearby massive E0 elliptical galaxies, NGC 1407, NGC 4486 (better known ad M87), and NGC 5846, belonging to the SLUGGS survey [70][71][72], a spectrophotometric survey of more than 4000 extragalactic GCs around 27 early-type galaxies at low redshift.This test was essential to verify whether the boost of the gravitational field could be determined by the gravitational permittivity alone, without any dependence on the refraction of the RG gravitational field lines.The results of this work are also summarised in [54].
Since NGC 1407, NGC 4486, and NGC 5846 have a minor-to-major axis ratio q of 0.95, 0.86, and 0.92, and consequently, an ellipticity ε = 1 − q of 0.05, 0.14, and 0.08, they can be approximated as spherical systems.Cesare et al. [47] chose these galaxies since the kinematic information of SLUGGS galaxies was measured up to their most external regions (≳10 effective radii), where Newtonian gravity does not work anymore unless DM is present, and the viability of RG could thus be tested.It was possible to measure these extended kinematic profiles thanks to GCs, which are present in these galaxy regions.In general, the measurement of the kinematic information in the external regions of elliptical galaxies can be uniquely performed thanks to the presence of kinematic tracers, such as X-ray emitting gas [73], planetary nebulae [74], and GCs [70], which settle in the galaxies outskirts.Specifically, the choice of [47] was driven by the fact that the GC population divides in two subpopulations of blue and red GCs, having different kinematic properties and formation histories.Specifically, the kinematics of the red GCs is usually similar to that of the stars in the host galaxy, maybe due to a similar formation history.Instead, the velocity dispersion of the blue GCs tend to be larger than that of the red GCs [70].In this way, it was possible to constrain RG from two, rather than one, extended kinematic profiles, providing a more robust test for the theory.

Mass Model
Cesare et al. [47] modelled, at the same time and for each E0 galaxy, the root-meansquare (RMS) velocity dispersions of the stars, concentrated within about one effective radius (R e ; the radius which encloses half the galaxy's luminosity), and of the blue and red GCs, distributed up to the galaxy outskirts.In particular, the kinematic profiles of the GCs were ∼5, ∼48, and ∼13 times more extended than the kinematic profiles of the stars in NGC 1407, NGC 4486, and NGC 5846.Cesare et al. [47] derived the model for the RMS velocity dispersion from the RG gravitational field computed with Equation (4) for spherical systems.To model the total mass M(< r) enclosed within the spherical radius r present in Equation ( 4), Cesare et al. [47] considered the contribution of (1) the stars, (2), the X-ray emitting gas, and (3) the central super massive black hole (SMBH).The total mass of GCs only contributed ∼1% to the total galaxy mass and it was thus neglected in the computation of the gravitational field.However, their surface and 3D number densities were modelled, since they were needed for the dynamical model (see Section 4.2).
The quantities entering the dynamical model were differently modelled for the three galaxies.The surface brightness of the stars of NGC 1407, measured in the B-band [75], was modelled with a Sérsic profile [58] of this kind, with parameters from [70] (Equation ( 7) in [47]): where (Equation ( 8) in [47]) is such that half of the total luminosity of the stars is enclosed within 1 R e .The circularised radius R is given for every galaxy by (Equation ( 6) in [47]): where x ′ and y ′ are coordinates oriented along the major and minor axes of the galaxies.Integrating Equation (24) with the Abel integral below yields the 3D luminosity density distribution for the stars of NGC 1407 (Equation ( 12) in [47]): and integrating Equation ( 27) in spherical coordinates yields its cumulative stellar luminosity profile (Equation ( 14) in [47]): For the surface brightness of the stars of NGC 4486 and NGC 5846, measured in the r-band [76], Cesare et al. [47] adopted a multi-Gaussian expansion (MGE) model [77] in a spherical configuration, with the parameters of [76] (Equation (10) in [47]): The 3D luminosity density distribution of the stars of NGC 4486 and NGC 5846 was given by (Equation ( 13) in [47]) and their cumulative luminosity stellar profile L * (< r) was obtained by Equation (28), where ν(r ′ ) is Equation (30).
The cumulative mass profile of the stars in each galaxy, M * (< r), was obtained by (Equation ( 28) in [47]): where L * (< r) is Equation ( 28), and Υ is the stellar mass-to-light ratio, one of the seven free parameters of the dynamical model.The 2D number density profiles of the blue and red GCs in NGC 1407 and NGC 5846 were modelled with a Sérsic profile [58] of this kind (Equation ( 15) in [47]): where Cesare et al. [47] adopted the parameters of [75] for NGC 1407 and estimated the parameters from the data of [70] with an MCMC for NGC 5846.For NGC 4486, they adopted a different parametrisation (Equation ( 16) in [47]): with the parameters of [78].The 3D number density profiles were obtained by integrating Equations ( 32) and (33) with the Abel integral (Equation ( 17) in [47]): Cesare et al. [47] modelled the mass density profile of the gas of NGC 1407 and NGC 5846 using the two-β functional form (Equations ( 18) and ( 19) in [47]): with the parameters of [79] for NGC 1407, and estimating the parameters from the data of [80] with an MCMC for NGC 5846.In Equation ( 35), µ = 0.6 and m H = 1.66054 × 10 −27 kg are the mean molecular weight and the atomic unit mass, respectively.For NGC 4486, they adopted this profile with the parameters of [81] (Equation ( 20) in [47]): The cumulative mass density profile of the gas, M gas (< r), was obtained by integrating in spherical coordinates ρ gas (r) from Equation (35) for NGC 1407 and NGC 5846, and Equation (36) for NGC 4486 (Equation ( 14) in [47]): The central SMBH was modelled as a point mass with values of M • = 4.5 +0.9 −0.4 × 10 9 M ⊙ , M • = 6.2 +0.4 −0.5 × 10 9 M ⊙ , and M • = 1.1 +0.1 −0.1 × 10 9 M ⊙ for NGC 1407, NGC 4486, and NGC 5846, obtained from [82].
The gravitational field dϕ/dr in Equation ( 38) is given by the RG field in a spherical configuration (Equation ( 4)), and thus, Equation (38) transforms into (Equation ( 25) in [47]): The cumulative mass profile M(< r) in Equation ( 40) is given by the sum of the cumulative mass profiles of the stars M * (< r) (Equation ( 31)), the gas M gas (< r) (Equation ( 37)), and the SMBH (end of Section 4.1): taken from Equation ( 26) in [47].This dynamical model has seven free parameters.Four parameters are common to the three tracers, the stellar mass-to-light ratio Υ, and the three RG parameters, ϵ 0 , Q, and ρ c , entering the gravitational permittivity ϵ(ρ) in Equation (40), which generate the RG gravitational potential to which every kinematic tracer is subject.Three parameters are specific to each tracer, and they are the orbital anisotropy parameters β * , β B , and β R .
To explore the parameter space, Cesare et al. [47] adopted a MCMC method based on a Metropolis-Hastings acceptance criterion.They assumed a uniform prior on every free parameter.For Υ, the uniform prior was based on the SPS models of Humphrey et al. [88] and Zhang et al. [79] for NGC 1407, and of Bell et al. [89] and Zibetti et al. [90] for NGC 4486 and NGC 5846.The space of the three RG free parameters was explored in the ranges of (0.00, 1.00], [0.01, 2.00], and [−27.00,−23.00], for ϵ 0 , Q, and log 10 [ρ c (g/cm 3 )], respectively.At last, they explored the parameter space of the parameters B * = − log 10 (1 in the uniform range of [−1.5, 1.0], which spans from very tangential to very radial orbits. The model given by Equation (40) was compared with the kinematic measurements of [75] for NGC 1407, and of the ATLAS 3D survey [91] for NGC 4486 and NGC 5846, to constrain the free parameters.

Results
The RG model (40) properly interpolated the kinematic profiles of the three tracers in each E0 galaxy with sensible parameters (Figure 6).Root-mean-square velocity dispersion profiles computed with RG with the parameters resulting from the MCMC analysis of [47], considering the stars, blue GCs, and red GCs' kinematic profiles at the same time for each elliptical E0 galaxy.(Top panels) the green, blue, and red solid lines represent the RG models for the stars, blue GCs, and red GCs, respectively.The dashed lines with the same colours are computed with the same Υ and β parameters as the solid lines but with Newtonian gravity without DM.(Bottom panels) zoom-in of the upper panels for the stars alone.Green solid and dashed lines have the same meaning as in the top panels.Both in the top and in the bottom panels, the black, blue, and red dots with error bars show the kinematic data for the stars, blue GCs, and red GCs, respectively.The figure is reproduced from Figure 8 in [47].Credit: Cesare V., Diaferio A., and Matsakos T., A&A, 657, A133, 2022, reproduced with permission © ESO.
The reasons for the ϵ 0 tension can be multiple.A first possible reason could be due to the approximate procedure with which the unique combination of RG parameters in the DMS was derived.A limited freedom was given to this unique combination of parameters to vary, since the mass parameters, Υ and h z , were constrained to the values derived from a previous analysis and this led to a very small width of the posterior distribution of ϵ 0 (see the light blue shaded area and yellow, red, and black contours in the bottom-left panel of Figure 4).As mentioned in Section 3.3, a more complete analysis could be performed by exploring the 63-dimensional parameter space of the 2 × 30 mass parameters and the three RG parameters in the DMS.
If, instead, the posterior distribution of ϵ 0 is not particularly underestimated, other reasons could explain the ϵ 0 discrepancy, such as a too approximated model for the elliptical galaxies, a wrong functional form for the gravitational permittivity ϵ(ρ), or a fundamental problem with RG.One of the approximations introduced in the analysis of elliptical galaxies was to consider these systems as isolated, whereas these kinds of galaxies typically live in dense environments.Specifically, NGC 1407 and NGC 5846 settle within galaxy groups and NGC 4486 is the central galaxy of the Virgo cluster.Differently from Newtonian gravity, the RG gravitational field of a system, as a MOND gravitational field, depends on the gravitational field produced by the environment where this system settles.Therefore, this environmental effect might explain, in part or entirely, the inconsistency between the ϵ 0 s.Another approximation adopted in modelling the elliptical galaxies was to neglect their net rotation, although these three galaxies are slow rotators [47].A caveat should be given since these elliptical galaxies were approximated as spherical systems due to their low minor-to-major axis ratio q.However, they could be significantly flattened oblate or prolate spheroidal systems seen along the symmetry axis.Moreover, or as an alternative, the functional form of the gravitational permittivity (Equation ( 10)) might not be suitable, in particular in the low-density regime, and a more precise answer in this sense might be given by the weak field limit (WFL) of the covariant formulation of RG [46].If all these reasons are not able to alleviate or cancel the ϵ 0 tension, this might be due to a fault in the theory that needs to be fixed.To conclude, a remarkable fact is that the best agreement between elliptical and disk galaxies was obtained for the ρ c parameter, which sets the density scale where the transition between the Newtonian and the RG regimes occurs, as a 0 in MOND.

Dynamics of Galaxy Clusters
The encouraging results obtained for the dynamics of disk and elliptical galaxies suggest that RG is able to properly describe the dynamics on a galaxy scale.Therefore, a class of tests at a larger scale need to be performed for a more complete investigation of the theory.Some preliminary studies have already been performed by Matsakos and Diaferio [45] to study the viability of RG in modelling the hot X-ray-emitting gas temperature radial profile of galaxy clusters.Matsakos and Diaferio [45] modelled the gas temperature radial profiles of two low-redshift and relaxed galaxy clusters, A1991 and A1795, with two different gas temperatures (Figure 7).Specifically, the two clusters have a spectroscopic gas temperature, averaged between 70 kpc and r 500 , of 2.83 keV (A1991) and 9.68 keV (A1795), where r 500 is the distance from the centre of the cluster where the average density is 500 times larger than the critical density of the Universe if DM is assumed to be present [92].These data are taken with the Chandra satellite [92].The model that Matsakos and Diaferio [45] adopted presented different assumptions.They assumed dynamical equilibrium at every distance from the cluster centre and a spherical configuration for the density distribution, deriving the equation of hydrostatic equilibrium from the RG gravitational field obtained, in turn, from Equation (4).They also assumed the validity of the ideal-gas equation of state and the smooth step function given by Equation ( 10) for the gravitational permittivity with RG parameters fixed a priori and not left free to vary.In particular, they set the RG parameters to ϵ 0 = 0.045 and 0.065 for A1991 and A1795, respectively, Q = 2, and ρ c = 10 −24 g cm −3 .At last, they assumed in A1991 the presence of a stellar component of mass 5 × 10 10 M ⊙ up to a distance from the cluster centre of R < 10 kpc.This model resulted in an agreement with the Chandra data of the two galaxy clusters.Figure 7 shows the RG models as black solid lines and the Newtonian expectations calculated with the same mass models as in RG, without considering any DM presence, as black dashed lines.This model is certainly simplistic and involves a too small sample of galaxy clusters.Therefore, it certainly needs to be extended to understand if RG is able to correctly describe the gas temperature radial profiles of galaxy clusters, and further tests have to be performed to investigate the capacity of RG to model the dynamics of these systems.An important test that has to be performed is to repeat the above analysis by constraining the RG parameters from the temperature data, and not by fixing them a priori, to verify if the obtained parameters are in agreement with the values obtained on a galaxy scale, to better test their universality.
Another aspect of RG on a cluster scale that Matsakos and Diaferio [45] began to explore is the gravitational interaction between galaxies in galaxy groups and clusters.These interactions are not simple to handle in RG, since the field lines in each galaxy suffer from refraction in low-density regions distant from their centres.A first-order analysis could be performed by assuming that the gravitational field at large distances from the centre of each galaxy decreases proportional to R −1 and that the only difference with respect to Newtonian gravitational field is that it is enhanced by ϵ −1 0 .With these assumptions, the BTFR (Equation ( 8)) can be extended to galaxy groups and clusters.In particular, the following expression for the BTFR would hold for galaxy groups and clusters (Equation (6.1) in [45]): where k is the number of galaxies, b is the BTFR normalisation, and f < 1 is a geometric factor that describes that, accounting for the anisotropic geometry of the RG field in disk galaxies, only a fraction of the k galaxies in the clusters will cross the plane of a specific disk galaxy where the field is ∝ R −1 .In Equation ( 42), v f represents the circular velocity of the galaxies in a specific galaxy group or cluster with total baryonic mass M (mainly made by the intracluster medium (ICM)), whereas in Equation ( 8), v f represents the circular velocity of the stars in a specific galaxy with total baryonic mass M. Figure 8 represents the BTFR extended to galaxy groups and clusters (Equation ( 42), black dashed line) for f /ϵ 0 = 1.3 and the BTFR of galaxies (Equation ( 8), black solid line).The black dots in the bottom-left part of the plot represent data points for galaxies taken from [93], having a mass range of M ∼ [10 9 , 10 12 ] M ⊙ and the black dots in the top-right part of the plot represent data points for galaxy groups and clusters taken from [94,95], having a mass range of M ∼ [10 12 , 10 15 ] M ⊙ .We can see that the two sets of data points distribute around the corresponding BTFR model.The open circles and squares represent simulated galaxies and they distribute around Equation (8) (black solid line).The observed normalisation of the BTFR extended to galaxy groups and clusters is smaller than the normalisation of the BTFR of galaxies but its slope is the same as for Equation (8), consistent with the approximate RG expectations.The comparison of this extended BTFR with the observed one for galaxy groups and clusters would represent another test for RG and a method to set a lower bound for ρ c and to constrain ϵ 0 .
Figure 8.The BTFR up to mass and velocity scales to include galaxy groups and clusters, besides galaxies.The black solid line is the BTFR given by Equation ( 8), and the black dashed line is the BTFR given by Equation (42) for f /ϵ 0 = 1.3.The open circles and squares represent simulated galaxies [45].The black dots that concentrate in the bottom-left part of the plot are observational data from [93] and refer to galaxies.The black dots that concentrate in the top-right part of the plot are observational data from [94,95] and refer to galaxy groups and clusters.The two set of measurement points concentrate around the corresponding BTFR model.The figure is reproduced from Figure 17

Covariant Refracted Gravity
A covariant version of RG was recently formulated by Sanna et al. [46].The fact that in RG, the transition between the Newtonian and the modified gravity regimes is regulated by a scalar quantity, the mass density ρ, rather than by a vector quantity, such as the acceleration as in MOND, made the building of a covariant extension of RG less challenging than in MOND (e.g., [96]).Some relativistic extensions of MOND, such as tensor vector scalar gravity (TeVeS) [96], were formulated while presenting some problems [97,98], even if more recent results might look more promising [99,100].
Covariant refracted gravity (CRG) is formulated as a scalar tensor theory with the presence of a single scalar field, φ, nonminimally coupled to the metric, which accounts for the phenomenologies of both DM on a galaxy scale and DE on larger scales, i.e., explaining the accelerated cosmic expansion.This peculiar feature of a unified dark sector is shared by a restricted class of modified theories of gravity (e.g., [101][102][103][104][105][106][107][108][109][110][111][112][113]) and is suggested by the a 0 ∼ H 0 ∼ √ Λ observed intriguing coincidence [46] (see Section 1).The general action of scalar-tensor theories is (Equation (3) in [46]): where the functional forms of the self-interaction potential V (φ) and the general differentiable function W (φ) of the scalar field φ define a specific scalar-tensor theory.For CRG, V (φ) and W (φ) are (Equation (7) in [46]): and (Equation ( 6) in [46]) where Ξ is a constant.Replacing Equations ( 44) and (45) in Equation ( 43), the CRG action becomes: In the WFL of the theory, the scalar field φ becomes twice the gravitational permittivity ϵ(ρ).Indeed, the WFL of CRG yields the original classical formulation of RG (Equation (24) in [46]): which is the RG Poisson equation (Equation ( 1)) if φ = 2ϵ(ρ).This result confirms that the scalar field mimics the phenomenology of DM on a galaxy scale.The fact that in RG, the modification of the law of gravity depends on a density rather than on an acceleration scale might be less intuitive, since the acceleration scale a 0 emerges from several pieces of evidence on a galaxy scale, such as the BTFR, the MDAR, and the RAR (see Section 1).Yet, the acceleration scale a 0 emerges from the WFL of CRG.This can be seen by calculating the CRG gravitational field far from a spherical source with density ρ s (r), monotonically decreasing with r, settling in a homogeneous background with constant density ρ bg .Sanna et al. [46] found that, far from this source, the transition from Newtonian to RG regimes occurs when the acceleration dϕ/dr goes below the acceleration scale (Equation ( 28) in [46]): where ρ(r) = ρ s (r) + ρ bg .This recalls the MOND theory, where the acceleration scale a 0 regulates the gravity behaviour.At large distances from the spherical source, the following limit holds (Section 3.4.1 in [46]): and therefore, Equation ( 48) reduces to (Section 3.4.1 in [46]): Since Ξ ∼ Λ, as found by independent calculations [46] explained below, Equation (50) implies that a Ξ ∼ 10 −10 m s −2 , which coincides with the MOND acceleration scale a 0 .
Sanna et al. [46] also derived with CRG the modified Friedmann equations for a flat, expanding, homogeneous, and isotropic Universe described with the Friedmann-Lema ître-Robertson-Walker (FLRW) metric, where the Universe content is modelled as a perfect fluid.In these modified Friedmann equations, the term (Equation ( 41) in [46]) appears, where H(t) is the Hubble parameter.This term is analogous to the density parameter related to the cosmological constant Λ in ΛCDM (Section 4.2 in [46]): which indicates that Ξ ∼ Λ, i.e, that Ξ plays the role of Λ in ΛCDM, accounting for the accelerated expansion of the Universe.At the present time, Equations ( 51) and ( 52) become: and where H 0 is the Hubble parameter at the present time.We have seen before that Ξ sets the value of the acceleration scale a Ξ (Equation ( 48)), which defines the transition between Newtonian and RG regimes far from a spherical source, playing the role of a 0 in MOND and thus accounting for DM phenomenology.Since Ξ ∼ Λ, Ξ also accounts for DE phenomenology, providing a unification of the two dark sectors.Moreover, inserting Ξ in Equation (48), the observed relation a 0 ∼ √ Λ naturally emerges in CRG.By rearranging the modified Friedmann equations in CRG, Sanna et al. [46] found two solutions for the Hubble parameter H(t), referred to as CRG− and CRG+, and properly integrating H(t), they found the corresponding two solutions for the scale factor a(t). Properly combining CRG− and CRG+, the bound Ξ > 0 was derived.From the derived a(t), the luminosity distance D L could be calculated and constrained from the observed Hubble diagram of SNe Ia at high redshift z of the Supernova Cosmology Project Union 2.1 Compilation [114] (open circles with error bars in Figure 9).Comparing the CRG model with the data, the cosmological parameters could be derived to verify if the tensions observed in ΛCDM could be either reduced or cancelled, which would provide a fundamental test for CRG.By comparing the distance modulus µ = m − M, where m and M are the apparent and absolute magnitudes, respectively, derived from the luminosity distance D L computed in CRG for the solutions CRG− and CRG+ (dashed and solid lines in Figure 9) with the (z, µ) data of the high-redshift SNe Ia from [114], Sanna et al. [46] found that Ω Ξ,0 = 0.650 +0.005 −0.650 for CRG−, and Ω Ξ,0 = 0.650 +0.005 −0.085 for CRG+, assuming the values H 0 = 67.7 km s −1 Mpc −1 and Ω 0 = 0.31 for the Hubble constant and the ratio between the densities of the baryonic matter and the critical density of the Universe at the present epoch, respectively [1].Replacing the values of Ω Ξ,0 = 0.65 and H 0 = 67.7 km s −1 Mpc −1 in Equation ( 54) implies again that Ξ ∼ Λ, as independently determined before.
Sanna et al. [46] derived the equation of state of the effective DE in CRG, w DE = p DE /ρ DE , where p DE and ρ DE are the pressure and the density of the effective DE, by properly rearranging the CRG's Friedmann equations and by comparing them with the Friedmann equations of a general scalar-tensor theory with a nonminimal coupling between the scalar field and the metric.The w DE parameter depends on redshift z.At the present time, z = 0, the w DE parameter becomes (Equation ( 73) in [46]): Inserting w DE = −1, consistent with the equation of state of DE in ΛCDM [1], and Ω 0 ∼ 0.3, Ω Ξ0 ∼ 0.64 is obtained, which is consistent with the CRG−solution.
Values of w DE different from −1 are consistent with several pieces of evidence anyway (e.g., [115][116][117][118][119][120]).Generally, the observational constraints on the equation of state of the effective DE can depend on the adopted model, even if model-independent reconstructions exist [120].For CRG, the parametrisation [121,122] (Section 4.3 in [46]): can be assumed for w DE .At the present time (z = 0), Equation ( 56) only depends on w 0 , allowing a wide range of DE models, either with w DE < −1 (phantom models) or with w DE ≥ −1 [115][116][117]123,124].The parameter w 0 is constrained to be approximately in the range of w 0 ∈ [−1.18, −0.85] by measurements of the baryonic acoustic oscillation (BAO), SNe Ia, and CMB [118,125].This limit on w 0 translates into a limit on Ω Ξ0 , in agreement with the values of Ω Ξ0 estimated from the SNe Ia data.

Discussion and Conclusions
In this work, the main analyses and results of the theory of modified gravity RG were summarised.RG demonstrated the ability to model the dynamics of disk and elliptical galaxies with sensible mass and structural parameters (stellar mass-to-light ratios, disk scale heights, and orbital anisotropy parameters) and with RG parameters consistent among the different galaxies, suggesting their universality.Preliminary encouraging results were obtained at the scale of galaxy clusters and a covariant extension of the theory might look promising, since it seems to properly describe the accelerated expansion of the Universe, to retrieve the MOND acceleration scale a 0 at a galaxy scale, and to explain both the DM and DE sectors with a single scalar field.
Several further studies could be performed to complete the tests of this gravity theory.On a galaxy scale, two issues presented by the theory require additional investigation.They are (1) the prediction of a RAR with a too large intrinsic scatter and correlations between its residuals from Equation ( 22) and some galaxy properties, and (2) the tension presented by the vacuum permittivity ϵ 0 , which might indicate its nonuniversality.Matsakos and Diaferio [45] showed that RG might be able to reproduce the BTFR and the MDAR.Instead, the RAR built with RG presents some strong correlations, at more than 5σ, with some radially dependent galaxy properties (the radius, the stellar surface density profile, and the gas fraction profile).However, the RAR built from DMS data also shows some, albeit weaker, correlations with some galaxy properties, and a further investigation is needed to understand whether the correlations presented by the RG RAR are partially driven by the data correlations.The correlations presented by the DMS data are apparently at odds with the claimed uncorrelations observed in the RAR derived from SPARC data [67], suggesting a difference between the two samples.Moreover, the question of the RAR is even more puzzling, since the DMS is not the only sample where a correlation between the RAR residuals and some galaxy properties has been observed.Di Paolo et al. [126] found a correlation between the RAR residuals and the galaxy radius from the accurate mass profiles of 36 dwarf disk spirals and 72 LSB galaxies.A better assessment of whether the result obtained in RG for the RAR depends on the theory itself or on the chosen galaxy sample can be performed by reproducing in RG the RAR of SPARC disk galaxies [67], which, differently from DMS galaxies, are nearly edge-on and thus, their measured rotation curves, from which the RAR is derived, are much more accurate.Moreover, the SPARC sample counts much more galaxies than the DMS (153 vs. 30 galaxies), which would make this study much more statistically significant.For this study, the rotation curves of SPARC galaxies have to be modelled with RG, and the RAR has to be built from these models.The scatter and the residuals from Equation ( 22) of the obtained RAR have to be compared with the results of Lelli et al. [69].Another interesting study to be performed with the RAR would be to compute the RAR in RG separately for the groups of normal spirals and dwarf galaxies present in SPARC and in the sample considered by Di Paolo et al. [126].This study might shed light on the works of Santos-Santos et al. [127] and Di Paolo et al. [126], which found that the RAR might differently behave in these two categories of galaxies, in addition to providing a further test for RG theory.
As said at the end of Section 4.3, the tension on the vacuum permittivity ϵ 0 might have different origins, either the approximate derivation of the unique ϵ 0 in the DMS sample, a too simplistic dynamical model for the elliptical galaxies, an incorrect assumption for the functional form of the gravitational permittivity ϵ(ρ) (Equation ( 10)), or, in the worst case, a fundamental issue with the RG theory.Generally, the study of elliptical galaxies is incomplete, both because the considered sample only counted three galaxies and because it only considered E0 ellipticals, with a nearly spherical shape.To better investigate the possible (non)universality of ϵ 0 , the performed study on the three E0 galaxies should be repeated removing the adopted approximations.In the new model, the interaction with the surrounding environment of the ellipticals should be considered, possibly resorting to N-body simulations.Additional E0 galaxies with extended kinematic profiles should be considered to enrich the galaxy sample in this new study.The E0 galaxies in the ePN.S survey [74], where extended kinematic profiles up to ∼13 R e from the galaxy centres were measured thanks to the presence of planetary nebulae in the galaxies outskirts, should be suitable candidates.Other candidates ideal to test the viability of RG are the round elliptical galaxies in the samples of [128][129][130].Their velocity dispersion profiles present a flattening beyond a certain radius, where both this radius and the intensity of the velocity dispersion in this plateau are in agreement with MOND expectations.In the outer regions of spherical systems, where the density drops below ρ c , RG predicts a field proportional to R −2 , as in the Newtonian case (see Section 2).Modelling in RG the velocity dispersion profiles of the round elliptical galaxies in [128][129][130] would show whether they can be reproduced by a ∝ R −2 field in the galaxy outskirts or if they necessarily require a ∝ R −1 field as in MOND, which would represent an issue with the current formulation of RG.
After focussing on spherical elliptical galaxies alone, a more complete study should be performed considering elliptical galaxies with different ellipticities, to further test the capability of RG of modelling the dynamics of systems with different shapes with a unique set of RG parameters.This study would also permit one to investigate the positive correlation found by Deur [49,50] between the total mass-to-light ratios and the ellipticities of elliptical galaxies, already mentioned in Section 1, which naturally arises within the RG context.Employing elliptical galaxies with kinematic profiles ∼10 times more extended than the data used by [49,50] might be crucial to validate or reject this correlation.All these studies might shed light on the correctness of Equation (10) for the gravitational permittivity ϵ(ρ).
A further test of the RG prediction that the flatter the system, the larger the mass discrepancy, if considered in Newtonian gravity, might be performed with the kinematic data of LSB galaxies (dwarf and dSph galaxies) and GCs, in addition to elliptical galaxies with different flatness degrees.Indeed, dwarf galaxies and GCs have similar baryonic masses but very different dynamical properties, the former being one of the darkest systems on a galaxy scale and the latter being nearly DM-free (Section 1).This last piece of evidence represented a challenge for MOND (e.g., [32][33][34][35][36][37][38]).This feature of RG might be tested by modelling the rotation curves of a sample of dwarf galaxies, e.g., in the Milky Way halo [131], in the LITTLE THINGS survey [132], in SPARC [67], and in Di Paolo et al. [126], and the internal velocity dispersion of a sample of GCs located in the most external regions of the Milky Way (e.g., [32][33][34][35][36][37][38]), where the background density drops below ρ c .If RG properly models the dynamics of these systems, it could explain why the dwarf galaxies are so "dark" if interpreted in Newtonian gravity, which can represent a challenge for ΛCDM [30] (See Section 1).Moreover, Matsakos and Diaferio [45] observed in a simplified framework that galaxies with a larger surface brightness in RG showed a smaller mass discrepancy, if interpreted in Newtonian gravity, than galaxies with a smaller surface brightness, which can naturally explain the difference between HSB and LSB galaxies and provide a solution for the cusp/core problem of ΛCDM (Section 1).An interesting study to test the viability of RG could be performed by modelling the dynamics of some dSph galaxies which present a high DM fraction, if interpreted in Newtonian gravity, but a not-so-squeezed shape, having an intrinsic ellipticity around 0.1-0.3[133,134].
On the scale of galaxy clusters, the tests are incomplete, and a deeper investigation is required.The sample considered by Matsakos and Diaferio [45] is too small, only counting two galaxy clusters, and the analysis of the temperature profiles is approximated.For their study, Matsakos and Diaferio [45] considered the data from [92], whose sample contains 11 other low-redshift and relaxed galaxy clusters.This analysis should be extended to these other galaxy clusters, by directly constraining the RG parameters from their temperature profiles to check their agreement with one another and with the results obtained on a galaxy scale.In a second step, the same analysis should be repeated by adopting a less approximate modelling by including the mass profiles of the individual galaxies present in the clusters, exploiting the results from N-body simulations.
In addition to the temperature profiles, RG should be tested on the dynamics of galaxy clusters.A study of this kind is already underway.Specifically, it should be verified whether RG is able to model the galaxy dispersion profiles of clusters of galaxies with sensible parameters.Possible samples for this analysis are CIRS [135] and HeCS [136], which present the radial velocity dispersion profiles taken from the Sloan Digital Sky Survey and the ROSAT All-Sky Survey of 130 galaxy clusters in a redshift range of [0.0, 0.3] and with different dynamical properties, from nonmerging and relaxed to merging and dynamically active.This study allows us to investigate how the RG parameters are affected by the environment, either relaxed or interacting, extending to larger scales the possible future study of elliptical galaxies.Also, the study of the BTFR extended to galaxy groups and clusters (Equation ( 42)) could be expanded, by considering how the uncertainties on the equation of state and on the entropy profiles of the gas and the possible local deviations from hydrostatic equilibrium due to X-ray gas flows might affect the values of the RG parameters.
Another important aspect of galaxy clusters that has not been investigated yet with RG and that should be studied to obtain a robust test for the viability of this theory, is the capability of RG to explain the hydrostatic mass bias, that is, the discrepancy between the masses of galaxy clusters measured with the gravitational lensing and the hot, X-ray-emitting gas methods [137].It would also be interesting to explore whether the dependence of the RG gravitational field on the shape of the considered system might explain gravitational lensing effects.
Also, CRG leaves room for additional investigations.The emergence of the MOND acceleration scale a 0 in the WFL of CRG should be further explored.Indeed, differently from a 0 , a Ξ (Equation ( 48)) is not a constant but depends on the mass density of the source, even if this dependence drops at large distances from the source, since ρ s ≪ ρ bg .Future studies are required to test if the connection between a Ξ and a 0 subsists for a generic case besides a source with a specific density configuration.
The parameter w DE = p DE /ρ DE of the equation of state of the effective DE found by Sanna et al. [46] in CRG appears to be dependent on the redshift, differently from ΛCDM where w DE is a constant equal to −1.However, w DE tends to −1 at the present epoch, in agreement with the observations.The w DE parameter in CRG could be constrained from the measurements of the expansion rate of the Universe performed by the Dark Energy Survey (DES) (https://www.darkenergysurvey.org,accessed on 4 April 2024) [138], which has observed thousands of supernovae since 31 August 2013.Moreover, the Euclid mission (https://sci.esa.int/web/euclid,accessed on 4 April 2024) [139], launched on 1 July 2023, is expected to measure the variation in the cosmic acceleration to an accuracy better than the 10% level, which would be crucial to disentangle CRG from ΛCDM and other DE models [140].Future observations of the large-scale structure of the Universe and the evolution of DE will be crucial to constrain the value of w 0 in Equation ( 56) and consequently, the value of Ξ (see Section 6).
The Lagrangian density of CRG, derived from the CRG action (Equation ( 46)), should be constrained from the power spectrum of the temperature anisotropies of the CMB and the power spectrum of the matter density perturbations (e.g., [141,142]), for example, by comparing the CRG predictions with the latest measurements from the Planck satellite [143].By estimating the cosmological parameters from different low-redshift galaxy surveys, such as the Kilo-Degree Survey (KIDS) (https://kids.strw.leidenuniv.nl,accessed on 4 April 2024) [144], Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) (https://www.cfht.hawaii.edu/Science/CFHLS/,accessed on 4 April 2024) [145], and DES, and from the CMB power spectrum, we can also assess whether RG can solve the tensions observed in the cosmological parameters in the ΛCDM model.
Another essential test for CRG would consist in analysing the evolution of the density perturbations, at least in the linear regime, and how the scalar field and its perturbations would impact on the large-scale structure formation, evolution, and distribution (e.g., [146][147][148][149]).A possible way to accomplish this task, would be to modify the publicly available RAMSES code [150], to run N-body simulations in the CRG framework.The results from these N-body simulations could be compared with the data from the DES survey, which probed the formation of structures with weak gravitational lensing and galaxy clusters, and the distribution of galaxies with the two-point correlation function.Further constraints could be given by the measurements taken with the Dark Energy Spectroscopic Instrument (DESI) (https://www.desi.lbl.gov,accessed on 4 April 2024) [151], which started to observe in 2019, with Euclid, and with upcoming experiments such as the Square Kilometer Array (SKA) (https://www.skao.int/,accessed on 4 April 2024) [152].
To conclude, a fundamental test to disentangle RG theory from GR and other modified theories of gravity can be provided by a modelisation of the gravitational waves, first detected on 14 September 2015 by the LIGO experiment [153].Indeed, some differences between GR and modified theories of gravity can be better stigmatised in the context of linearised gravitation in this era of the recently started gravitational wave astronomy, as stressed and anticipated by [154].

Figure 1 .
Figure 1.Comparison between the behaviour of Newtonian (top panels) and RG (bottom panels) gravitational field lines for flat (left panels) and spherical (right panels) systems.The formulae reported for flattened systems, in the left part, are equivalent formulae for spherical systems to approximately illustrate the dependence of the field on the distance from the galaxy centre.The figure is reproduced from Figure 16 in [45].

Figure 3 .
Figure3.RG models (blue solid lines) of the rotation curves and the vertical velocity dispersion profiles against the DMS data (red dots with error bars) for three DMS galaxies (UGC 3091, UGC 3701, and UGC 4256).Panels (a): Rotation curves computed with the parameters found by modelling the rotation curves alone of each galaxy.Panels (b): Rotation curves computed with the parameters found by modelling the rotation curves and vertical velocity dispersion profiles of each galaxy at the same time.Panels (c): Vertical velocity dispersion profiles computed with the parameters found by modelling the rotation curves and vertical velocity dispersion profiles of each galaxy at the same time.The cyan solid lines are the RG vertical velocity dispersion profiles calculated with the same parameters as the blue curves but with disk-scale heights h z = h z,SR (Equation (21)), and the dashed magenta and green vertical lines show the bulge effective radius (R e in Equation (15)) and the bulge radius adopted in the disk-bulge decomposition for the surface brightness modelling, respectively.The figure is readapted from Figures D.1 and D.2 in[48].Credit: Cesare V., Diaferio A., Matsakos T., and Angus G., A&A, 637, A70, 2020, reproduced with permission © ESO.

Figure 4 .
Figure 4. Free parameters of the RG gravitational permittivity ϵ(ρ) estimated from the DMS and the three elliptical E0 galaxies.The purple squares with error bars represent the means of the permittivity parameters estimated at the same time from the rotation curves and the vertical velocity dispersion profiles of each DMS galaxy.The light blue shaded areas indicate the posterior distributions of the unique combination of RG parameters estimated from the entire DMS sample at the same time, where their medians and 1σ, 2σ, and 3σ levels are shown as green dots and yellow, red, and black contours, respectively.The orange squares with error bars are the permittivity parameters derived from each elliptical E0 galaxy.The figure is reproduced from Figure 9 in [47].Credit: Cesare V., Diaferio A., and Matsakos T., A&A, 657, A133, 2022, reproduced with permission © ESO.

Figure 5 .
Figure5.RAR models calculated with RG (blue solid lines) with the parameters obtained from the MCMC analysis of the rotation curves and vertical velocity dispersion profiles of each galaxy at the same time, superimposed onto the RAR of DMS data (red points with error bars).The black solid line is the RAR fitting relation (Equation (22)) obtained by McGaugh et al.[41] from SPARC galaxies and the black dashed line is the g obs = g bar relation, plotted as a reference.The figure is readapted from Figure11in[48].Credit: Cesare V., Diaferio A., Matsakos T., and Angus G., A&A, 637, A70, 2020, reproduced with permission © ESO.

Figure 6 .
Figure6.Root-mean-square velocity dispersion profiles computed with RG with the parameters resulting from the MCMC analysis of[47], considering the stars, blue GCs, and red GCs' kinematic profiles at the same time for each elliptical E0 galaxy.(Top panels) the green, blue, and red solid lines represent the RG models for the stars, blue GCs, and red GCs, respectively.The dashed lines with the same colours are computed with the same Υ and β parameters as the solid lines but with Newtonian gravity without DM.(Bottom panels) zoom-in of the upper panels for the stars alone.Green solid and dashed lines have the same meaning as in the top panels.Both in the top and in the bottom panels, the black, blue, and red dots with error bars show the kinematic data for the stars, blue GCs, and red GCs, respectively.The figure is reproduced from Figure8in[47].Credit: Cesare V., Diaferio A., and Matsakos T., A&A, 657, A133, 2022, reproduced with permission © ESO.

Figure 7 .
Figure 7. Emission-weighted projected temperature radial profiles of the hot, X-ray-emitting gas in the low-redshift and relaxed galaxy clusters A1991 (top panel) and A1795 (bottom panel).The solid and dashed lines are the RG and Newtonian models, respectively, and the dots with error bars are the measurements taken with the Chandra satellite from [92].The Newtonian expectations were calculated with the same mass model as in RG without considering the presence of DM.A Hubble constant of H 0 = 71 km s −1 Mpc −1 was adopted.The figure is reproduced from Figure 14 in [45].

Figure 9 .
Figure 9. Hubble diagram of SNe Ia modelled with CRG.The dashed and solid lines are the CRG− and CRG+ models, respectively.For the CRG− and CRG+ curves, the H 0 = 67.7 km s −1 Mpc −1 , Ω 0 = 0.31, and Ω Ξ,0 = 0.65 parameters were adopted.The open circles with error bars are the data from the Supernova Cosmology Project Union 2.1 Compilation [114].The figure is reproduced from Figure E.2 in [46].