Refracted Gravity Solutions from Small to Large Scales

1. Introduction

In the Universe, several pieces of evidence cannot reconcile with General Relativity (GR), if we only assume the presence of visible (baryonic) matter. A mass discrepancy of ∼80-90% is observed from the largest scales (e.g., the cosmic microwave background (CMB) [1], the large-scale structure [2,3], gravitational lensing effects where Bullet Cluster example is worth mentioning [4,5,6], light elements abundances [7], and galaxy clusters dynamics [8]) to the smallest scales (e.g., the flat trend of disk galaxies rotation curves [9,10,11]).

The most investigated solution to explain this phenomenology is to assume the presence of a non-baryonic and cold, i.e., non-relativistic 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 $\mathrm{\Lambda }$CDM, which assumes GR as 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 (SNaeIa) [13,14], and it can be identified with the cosmological constant $\mathrm{\Lambda }$.

Even if $\mathrm{\Lambda }$CDM can account for the majority of the observations in the Universe, it presents some problems, both on large and on 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 the early Universe [1,18,19]. On small scales, we observe several discrepancies between $\mathrm{\Lambda }$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).

It is remarkable to cite also the presence of some observed regularities on galaxy scale which are hard to explain in a $\mathrm{\Lambda }$CDM paradigm, where a stochastic merging process of structures is invoked. Some of these regularities are very tight scaling relations between a property of the dark matter and of the baryonic matter in galaxies, which might be counter-intuitive since DM represents the ∼90% of the galaxy content, whereas baryonic matter only represents its ∼10%. Among these relations, we can mention (1) the Baryonic Tully-Fisher Relation [22], (2) the Mass Discrepancy-Acceleration Relation [23], and (3) the Radial Acceleration Relation [24]. The three relations see the emergence of the same acceleration scale $a_0\simeq 1.2\times {10}^{-10}$ m s⁻², which, if expressed in natural units, is $a_0\sim H_0\sim \sqrt{\mathrm{\Lambda }}$, suggesting a connection between the DM and the DE sectors. These coincidences are even less intuitive in a $\mathrm{\Lambda }$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 ways is MOdified Newtonian Dynamics (MOND) [25,26,27]. 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 predicted 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 [28] 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 $\rho$, rather than of the acceleration a. The Poisson equation of RG is modified at first member by the presence of the gravitational permittivity, $ϵ\left(\rho \right)$, a monotonic increasing function of $\rho$ 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 [29] and it seems to describe both the DE and DM sectors with a single scalar field and to reproduce the the Hubble diagram of SNaeIa.

RG has obtained some encouraging results in modelling the dynamics of galaxies and galaxy clusters, and on the covariant side. The paper develops as follows. Section 2 describes RG theory. Section 3, Section 4, Section 5, and Section 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.

2. Refracted Gravity

Refracted Gravity (RG) is a theory of modified gravity formulated in a classical way, which can be interpreted in analogy to electrodynamics in matter [28]: 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:

$$\nabla ·\left[ϵ\left(\rho \right)\nabla \varphi \right]=4\pi G\rho ,$$

(1)

where $\varphi$ is the gravitational potential and $ϵ\left(\rho \right)$ is the gravitational permittivity. The following asymptotic limits are adopted for the gravitational permittivity:

$$ϵ\left(\rho \right)=\left\{\begin{array}{cc}1,\hfill & \rho \gg {\rho }_{\mathrm{c}}\hfill \\ {ϵ}_0,\hfill & \rho \ll {\rho }_{\mathrm{c}},\hfill \end{array}\right.$$

(2)

where $0<{ϵ}_0\le 1$ and ${\rho }_{\mathrm{c}}$ are the gravitational permittivity in vacuum and the critical density, respectively, two of the three free parameters of the theory. When $\rho \gg {\rho }_{\mathrm{c}}$, Newtonian Poisson equation is recovered:

$${\nabla }^2{\varphi }_{\mathrm{N}}=4\pi G\rho ,$$

(3)

whereas, when $\rho \ll {\rho }_{\mathrm{c}}$, we are in fully RG regime.

Differently from MOND, the RG gravitational field has a different behaviour for spherical and flat systems. For spherical systems, the RG Poisson equation (Equation (1)) reduces to:

$$\frac{\partial \varphi }{\partial r}=\frac{1}{ϵ\left(\rho \right)}\frac{GM(<r)}{r^2}=\frac{1}{ϵ\left(\rho \right)}\frac{\partial {\varphi }_{\mathrm{N}}}{\partial r},$$

(4)

where $M(<r)$ is the mass of the spherical system enclosed within the spherical radius r and $\partial {\varphi }_{\mathrm{N}}/\partial 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 $\rho \ll {\rho }_{\mathrm{c}}$, i.e., where $ϵ\left(\rho \right)\to {ϵ}_0$ (right panels of Figure 1).

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:

$$\frac{\partial ϵ}{\partial \rho }\nabla \rho ·\nabla \varphi +ϵ\left(\rho \right){\nabla }^2\varphi =4\pi G\rho .$$

(5)

In this configuration, the RG field depends both on the density $\rho$ (second term in the left-hand side of Equation (5)), as for spherical systems, and on its gradient (first term in the left-hand side of Equation (5)). In particular, the term “$\frac{\partial ϵ}{\partial \rho }\nabla \rho ·\nabla \varphi$” is different from zero in non-spherical 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, in RG is explained 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 (5), 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 [30,31] and the different DM quantity of globular clusters and dwarf galaxies, where the former are nearly spherical and DM-free (e.g., [32]) and the latter are flatter and one of the darkest objects in the Universe (e.g., [33]).

Where $\rho \ll {\rho }_{\mathrm{c}}$ in flat systems, e.g., in the outskirts of disk galaxies, the gravitational field assumes the asymptotic trend $\partial \varphi /\partial R\sim \sqrt{|\partial {\varphi }_{\mathrm{N}}/\partial R|a_0}\propto R^{-1}$, where $a_0$ is the MOND critical acceleration set by the normalization of the BTFR:

$$V_{\mathrm{f}}^4=G{a}_{0}M,$$

(6)

where $V_{\mathrm{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. In this regime, the RG field trend deviates from the Newtonian inverse square law. This asymptotic limit coincides with MOND asymptotic limit for the gravitational field where $a\ll a_0$, which might indicate that RG shares the majority of MOND successes on galaxy scale. The difference between RG and MOND is better observed in spherical systems, where the field trend in RG is Newtonian ($\propto R^{-2}$) where $\rho \ll {\rho }_{\mathrm{c}}$, preserving the Gauss theorem, but the field trend in MOND is $\propto R^{-1}$ where $a\ll a_0$, as in flattened systems.

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 $\rho$ was employed:

$$ϵ\left(\rho \right)={ϵ}_0+(1-{ϵ}_0)\frac{1}{2}\left\{tanh\left[\ln {\left(\frac{\rho }{{\rho }_{\mathrm{c}}}\right)}^Q\right]+1\right\},$$

(7)

represented in Figure 2. Besides ${ϵ}_0$ and ${\rho }_{\mathrm{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 (7) 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.

3. Dynamics of disk galaxies

Matsakos and Diaferio [28] presented some preliminary encouraging results, according to which RG properly models the BTFR of galaxies (Equation (6)) and the rotation curves of the high-surface brigthness (HSB) disk galaxy NGC6946 and low-surface brightness (LSB) disk galaxy NGC1560. Cesare et al. [34] demonstrated that RG can properly model the dynamics of disk galaxies. For their analysis, they considered the rotation curves and the vertical velocity dispersion radial profiles of 30 disk galaxies in the DiskMass Survey (DMS) [35] at low redshift. They chose this sample since these galaxies are close to face-on, which allowed the measurement of both their rotation curves and vertical velocity dispersions. In this way, Cesare et al. [34] 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 [36].

Cesare et al. [34] modelled with RG (1) the rotation curves alone of each DMS galaxy, (2) the rotation curves and the vertical velocity dispersions at the same time of each DMS galaxy, (3) the rotation curves and the vertical velocity dispersions of all the 30 DMS galaxies at the same time, and (4) the RAR of DMS galaxies.

3.1. Mass model

Cesare et al. [34] modelled the baryonic mass density profile $\rho$, which generates the total gravitational potential $\varphi$ since DM is not present, with (1) an axisymmetric exponential disk for the stars, (2) a spherical stellar bulge, and (3) two razor-thin disks for the atomic and molecular gas. The 3D mass density of the stellar disk is modelled with an exponential profile of this kind:

$${\rho }_{\mathrm{d}}(R,z)=\frac{\mathrm{{\rm Y}}{I}_{\mathrm{d}}\left(R\right)}{2h_z}exp\left(-\frac{\left|z\right|}{h_z}\right),$$

(8)

where $I_d\left(R\right)$ 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 $\mathrm{{\rm Y}}$ and $h_z$ are the stellar mass-to-light ratio and the disk scale height, two of the five free parameters of the dynamic model adopted by [34].

The surface brightness of the disk is 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 the “Renzo’s rule” [37]. The disk surface brightness in the innermost galaxy regions, where the bulge is dominant, is modelled with an exponential profile [38,39],

$$I_{\mathrm{d}}\left(R\right)=I_{\mathrm{d},0}exp\left(-\frac{R}{h_R}\right),$$

(9)

obtained by fitting the surface brightness points of the outermost stellar disk, where the stellar disk alone dominates.

DMS galaxies are disk-dominated, being the stellar disk-over-stellar bulge total luminosity ratio around $0.09$. For this reason, the bulge was approximated as a sphere, introducing negligible systematic errors. The 3D mass density of the bulge is modelled with this Abel integral:

$${\rho }_{\mathrm{b}}\left(r\right)=-\frac{\mathrm{{\rm Y}}}{\pi }{\int }_r^{+\infty }\frac{\mathrm{d}{I}_{\mathrm{b}}\left(R\right)}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}},$$

(10)

where $I_b\left(R\right)$ is the surface brightness radial profile of the bulge, and

$$r=\sqrt{R^2+z^2}$$

(11)

is the 3D radius. Being the bulge contribution subdominant, Cesare, et al. [34] assumed the mass-to-light ratio of the bulge equal to the mass-to-light ratio of the disk, $\mathrm{{\rm Y}}$, without introducing an additional degree of freedom. They modelled the surface brightness of the bulge with a Sérsic spherical profile [40]:

$$I_{\mathrm{b}}\left(R\right)=I_{\mathrm{e}}exp\left\{-7.67\left[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{1/n_s}-1\right]\right\}$$

(12)

To account for the seeing, which is not negligible since the measurements where taken from the ground with the 3.5 m diameter Calar Alto Observatory, they modelled the observed surface brightness of the bulge with Equation (12) convolved with a Gaussian point spread function.

The 3D mass density of the atomic and molecular gas are modelled with two razor-thin disks:

$${\rho }_{\mathrm{atom}}(R,z)={\Sigma }_{\mathrm{atom}}\left(R\right)\delta \left(z\right),$$

(13)

and

$${\rho }_{\mathrm{mol}}(R,z)={\Sigma }_{\mathrm{mol}}\left(R\right)\delta \left(z\right),$$

(14)

where ${\Sigma }_{\mathrm{atom}}\left(R\right)$ and ${\Sigma }_{\mathrm{mol}}\left(R\right)$ are the surface mass densities of the atomic and molecular gas, respectively, which Cesare, et al. [34] modelled with a linear interpolation of the data for the same reason as the stellar disk, and $\delta \left(z\right)$ is the Dirac $\delta$ function.

The total 3D mass density profile is given by summing Equations (8), (10), (13), and (14):

$$\rho (R,z)={\rho }_{\mathrm{d}}(R,z)+{\rho }_{\mathrm{b}}\left(r\right)+{\rho }_{\mathrm{atom}}(R,z)+{\rho }_{\mathrm{mol}}(R,z).$$

(15)

3.2. Dynamical model

From the total 3D mass density profile given by Equation (15), Cesare et al. [34] 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-Over-Relaxation method, given by combining the Jacobi and Gauss-Seidel methods [41]. From the RG potential, they derived the models for the rotation curve, $v\left(R\right)$, and the vertical velocity dispersion radial profile, ${\sigma }_z\left(R\right)$:

$$v(R,z=0)=\sqrt{R\frac{\partial \varphi (R,z)}{\partial R}},$$

(16)

and

$${\sigma }_z\left(R\right)=\sqrt{\frac{1}{h_z}{\int }_0^{+\infty }\left[{\int }_z^{+\infty }exp\left(-\frac{|z^{\prime }|}{h_z}\right)\frac{\partial \varphi (R,z^{\prime })}{\partial z^{\prime }}d{z}^{\prime }\right]dz},$$

(17)

where Equation (17) is the vertical velocity dispersion for a vertically-decreasing exponential disk whose density follows Equation (8).

The dynamical model has five free parameters: the stellar mass-to-light ratio $\mathrm{{\rm Y}}$, the stellar disk scale height $h_z$, and the three RG parameters, ${ϵ}_0$, Q, and ${\rho }_{\mathrm{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 Metropolis-Hastings acceptance criterion. They adopted the following priors for the free parameters of the model. For $\mathrm{{\rm Y}}$, they adopted a Gaussian prior based on the stellar population synthesis (SPS) model of Bell and de Jong [42], keeping into account that the surface brightness is measured in the K-band, where the Gaussian tail for $\mathrm{{\rm Y}}<0$ was set to 0. For $h_z$, they adopted a Gaussian prior peaked on $h_{z,\mathrm{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. [43]:

$${log}_{10}\left(\frac{h_R}{h_{z,\mathrm{SR}}}\right)=0.367{log}_{10}\left(\frac{h_R}{\mathrm{kpc}}\right)+0.708\pm 0.095,$$

(18)

where the term $\pm 0.095$ indicates the ∼25% intrinsic scatter. Estimating $h_{z,\mathrm{SR}}$ was essential to have a comparison reference for DMS galaxies where, due to their nearly face-on configuration, their $h_z$ cannot 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}\left[{\rho }_{\mathrm{c}}(\mathrm{g}/{\mathrm{cm}}^3)\right]$.

3.3. Results

RG models 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) are computed with the parameters estimated from the rotation curves of each galaxy alone. The blue curves for the rotation curves in panels (b) are computed with the parameters found by modelling the rotation curves and the vertical velocity dispersion profiles of each galaxy at the same time. 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 equal to $h_{z,\mathrm{SR}}$ (Equation (18)). In all panels, the dashed magenta vertical lines identify the bulge effective radius ($R_{\mathrm{e}}$ in Equation (12)), 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 models both kinematic profiles and that the substructures of the rotation curves are properly reproduced following the 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. [34] found mass-to-light ratios $\mathrm{{\rm Y}}$ consistent with the SPS model of Bell and de Jong [42]. In the only-rotation curve analysis, the estimated $h_z$ were generally consistent with the $h_{z,\mathrm{SR}}$ derived from Equation (18). Instead, in the combined analysis, the estimated $h_z$ systematically resulted smaller than the correspondent $h_{z,\mathrm{SR}}$. Angus, et al. [44] found a similar result by modelling the DMS galaxies with QUMOND theory [45], a modified gravity version of MOND theory. Yet, this result seems not to 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. [34] and by Angus, et al. [44] 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 [46]. The young stellar population distributes in a disk thinner and with a smaller vertical velocity dispersion ${\sigma }_z$ than the old stellar population [46]. Instead, the disk scale heights of edge-on galaxies, from which Equation (18) is derived to determine $h_{z,\mathrm{SR}}$, are directly measured from near-Infrared photometry, whose signal is dominated by the old stellar population [46], having higher $h_z$ and ${\sigma }_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. [34] computed this factor which resulted equal to $1.63$, in agreement with the work of Aniyan, et al. [46]. Cesare et al. [34] artificially increased of this factor the vertical velocity dispersions of five DMS galaxies, obtaining disk scale heights in agreement with the correspondent $h_{z,\mathrm{SR}}$.

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 distribution. This suggests their universality, as should be from the theory formulation. To better test their universality, Cesare, et al. [34] estimated a unique combination of RG free parameters from the two kinematic profiles of all the 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 are ${ϵ}_{0,\mathrm{Mean},\mathrm{DMS}}=0.56\pm 0.16$, $Q_{\mathrm{Mean},\mathrm{DMS}}=0.92\pm 0.71$, ${log}_{10}{\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]}_{\mathrm{Mean},\mathrm{DMS}}=-25.30\pm 1.22$ (purple squares with error bars in Figure 4), whereas the unique combination of free parameters is ${ϵ}_{0,\mathrm{Unique},\mathrm{DMS}}=0.{661}_{-0.007}^{+0.007}$, $Q_{\mathrm{Unique},\mathrm{DMS}}=1.{79}_{-0.26}^{+0.14}$, ${log}_{10}{\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]}_{\mathrm{Unique},\mathrm{DMS}}=-24.{54}_{-0.07}^{+0.08}$, where their estimates and error bars are 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$\sigma$, 2$\sigma$, and 3$\sigma$ levels. The code used to find this unique combination of RG parameters is written in C++ and parallelized with OpenMP. It is publicly available on GitHub under the name “astroMP”1 and it is described in [47,48].

¹https://github.com/alpha-unito/astroMP

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 provides a limited worsening of the agreement between the RG models and the kinematic data and that the unique ${ϵ}_0$, Q, and ${log}_{10}\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]$ agree with the correspondent RG parameters estimated from the single galaxies within $0.63$$\sigma$, $1.19$$\sigma$, and $0.62$$\sigma$, respectively, further enforces the hypothesis of the universality of these parameters. With a proper exploration of the 63-dimensional parameter space for this galaxy sample (one $\mathrm{{\rm Y}}$ 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. [34] also modelled the RAR of the DMS galaxies. McGaugh, et al. [24] fitted the RAR of 153 edge-on disk galaxies in the SPARC sample [50] with this relation:

$$g_{\mathrm{obs}}\left(R\right)=\frac{g_{\mathrm{bar}}\left(R\right)}{1-exp\left(-\sqrt{\frac{g_{\mathrm{bar}}\left(R\right)}{g_{†}}}\right)},$$

(19)

where $g_{\mathrm{obs}}\left(R\right)$ is the observed centripetal acceleration derived from the measured rotation curve $v_{\mathrm{obs}}\left(R\right)$,

$$g_{\mathrm{obs}}\left(R\right)=\frac{v_{\mathrm{obs}}{\left(R\right)}^2}{R},$$

(20)

$g_{\mathrm{bar}}\left(R\right)$ 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\pm 0.02\left(\mathrm{random}\right)\pm 0.24\left(\mathrm{systematic}\right)]\times {10}^{-10}$ m s⁻², which agrees with the MOND acceleration scale $a_0$ within 1$\sigma$.

The RG models of the RAR for the DMS galaxies (blue solid lines in Figure 5) properly interpolate the measured RAR of the DMS data (red dots with error bars in Figure 5) and reproduce the asymptotic limits of Equation (19). However, RG reproduces the RAR of DMS galaxies with a too large intrinsic scatter ($0.11$ dex vs $0.057$ dex as found by Li et al. [51] from SPARC galaxies) and with residuals strongly correlated, at largely more than 5$\sigma$, with some galaxy properties dependent on the distance from the galaxies centre, apparently at odds with the results found by Lelli et al. [52] with SPARC galaxies. Yet, this problem requires further investigation since also the RAR of DMS galaxies presents some, although weaker, correlations between its residuals from Equation (19) and some galaxy properties. A galaxy sample such as the DMS, made of only 30 galaxies and where the galaxies are close to face-on and, thus, rotation curves might not be so accurately measured, 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.

4. 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 systems, Cesare, et al. [49] 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 [53,54,55], 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 can 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 [36].

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 $\epsilon =1-q$ of $0.05$, $0.14$, and $0.08$, they can be approximated as spherical systems. Cesare, et al. [49] chose these galaxies since the kinematic information of SLUGGS galaxies was measured up to their most external regions (≳10 effective radii), where Newonian 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 done thanks to the presence of kinematic tracers, such as X-ray emitting gas [56], planetary nebulae [57] and GCs [53], which settle in the galaxies outskirts. Specifically, the choice of [49] was driven by the fact that the GC population divides in two sub-populations of blue and red GCs, having different kinematic properties and formation histories. In this way, it was possible to constrain RG from two rather than one extended kinematic profile, providing a more robust test for the theory.

4.1. Mass model

Cesare et al. [49] modelled, at the same time, the root-mean-square (RMS) velocity dispersions of the stars, concentrated within about one effective radius ($R_{\mathrm{e}}$; the radius which encloses half galaxy luminosity), and of the blue and red GCs, distributed up to the galaxy outskirts, of each E0 galaxy. In particular, the kinematic profiles of the GCs are ∼5, ∼48, and ∼13 times more extended than the kinematic profiles of the stars in NGC 1407, NGC 4486, and NGC 5846, respectively. Cesare et al. [49] 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. [49] 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 [58], was modelled with a Sérsic profile of this kind, with parameters from [53]:

$$I_{*}\left(R\right)=I_{\mathrm{e}}exp\left\{-b_{n_{\mathrm{s}}}\left[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{1/n_{\mathrm{s}}}-1\right]\right\}$$

(21)

where

$$b_{n_{\mathrm{s}}}\approx 2n_{\mathrm{s}}-\frac{1}{3}+\frac{4}{405n_{\mathrm{s}}}+\frac{46}{25515n_{\mathrm{s}}^2}$$

(22)

is such that half of the total luminosity of the stars is enclosed within 1 $R_{\mathrm{e}}$. The circularised radius R is given for every galaxy by:

$$R=\sqrt{q{x}^{\prime 2}+\frac{y^{\prime 2}}{q}},$$

(23)

where $x^{\prime }$ and $y^{\prime }$ are coordinates oriented along the major and minor axes of the galaxies.

Integrating Equation (21) with the Abel integral below yields the 3D luminosity density distribution for the stars of NGC 1407:

$${\nu }_{*}\left(r\right)=-\frac{1}{\pi }{\int }_r^{+\infty }\frac{\mathrm{d}{I}_{*}\left(R\right)}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}},$$

(24)

and integrating Equation (24) in spherical coordinates yields its cumulative stellar luminosity profile:

$$L_{*}(<r)=4\pi {\int }_0^r{\nu }_{*}\left(r^{\prime }\right)r^{\prime 2}\mathrm{d}{r}^{\prime }.$$

(25)

For the surface brightness of the stars of NGC 4486 and NGC 5846, measured in the r-band [59], Cesare, et al. [49] adopted a Multi-Gaussian Expansion (MGE) model in a spherical configuration, with the parameters of [59]:

$$I_{*}\left(R\right)=\sum _{k=1}^N\frac{L_k}{2\pi {\sigma }_k^2}exp\left(-\frac{R^2}{2{\sigma }_k^2}\right).$$

(26)

The 3D luminosity density distribution of the stars of NGC 4486 and NGC 5846 is given by

$${\nu }_{*}\left(r\right)=\sum _{k=1}^N\frac{L_k}{{\left(\sqrt{2\pi }{\sigma }_k\right)}^3}exp\left(-\frac{r^2}{2{\sigma }_k^2}\right)$$

(27)

and their cumulative luminosity stellar profile $L_{*}(<r)$ is obtained by Equation (25), where $\nu \left(r^{\prime }\right)$ is Equation (27). The cumulative mass profile of the stars, $M_{*}(<r)$, is obtained by:

$$M_{*}(<r)=\mathrm{{\rm Y}}{L}_{*}(<r),$$

(28)

where $L_{*}(<r)$ is Equation (25) and $\mathrm{{\rm Y}}$ 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 are modelled with a Sérsic profile of this kind:

$$N_{\mathrm{GC}}\left(R\right)=N_{\mathrm{e}}exp\left\{-b_{n_s}\left[{\left(\frac{R}{R_{\mathrm{e}}}\right)}^{1/n_s}-1\right]\right\},$$

(29)

where Cesare, et al. [49] adopted the parameters of [58] for NGC 1407 and estimated the parameters from the data of [53] with a MCMC, for NGC 5846. For NGC 4486, they adopted a different parametrisation:

$$N_{\mathrm{GC}}\left(R\right)=N_{0}exp\left[-{\left(\frac{R}{R_{\mathrm{s}}}\right)}^{\frac{1}{m}}\right],$$

(30)

with the parameters of [60]. The 3D number density profiles are obtained by integrating Equations (29) and (30) with the Abel integral:

$${\nu }_{\mathrm{GC}}\left(r\right)=-\frac{1}{\pi }{\int }_r^{+\infty }\frac{\mathrm{d}{N}_{\mathrm{GC}}\left(R\right)}{\mathrm{d}R}\frac{\mathrm{d}R}{\sqrt{R^2-r^2}}.$$

(31)

Cesare, et al. [49] modelled the mass density profile of the gas of NGC 1407 and NGC 5846 using the two-$\beta$ functional form:

$${\rho }_{\mathrm{gas}}\left(r\right)=\mu m_{\mathrm{H}}\sqrt{n_{\mathrm{g},1}^2{\left[1+{\left(\frac{r}{R_{\mathrm{c},1}}\right)}^2\right]}^{-3{\beta }_1}+n_{\mathrm{g},2}^2{\left[1+{\left(\frac{r}{R_{\mathrm{c},2}}\right)}^2\right]}^{-3{\beta }_2}},$$

(32)

with the parameters of [61], for NGC 1407, and estimating the parameters from the data of [62] with a MCMC for NGC 5846. In Equation (32), $\mu =0.6$ and $m_{\mathrm{H}}=1.66054\times {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 [63]:

$${\rho }_{\mathrm{gas}}\left(r\right)=\frac{{\rho }_0}{{(1+b^{\prime }{r}^2+c^{\prime }{r}^4+d^{\prime }{r}^6)}^{n^{\prime }}}.$$

(33)

The cumulative mass density profile of the gas, $M_{\mathrm{gas}}(<r)$, is obtained by integrating in spherical coordinates ${\rho }_{\mathrm{gas}}\left(r\right)$ from Equation (32), for NGC 1407 and NGC 5846, and Equation (33), for NGC 4486:

$$M_{\mathrm{gas}}(<r)=4\pi {\int }_0^r{\rho }_{\mathrm{gas}}\left(r^{\prime }\right)r^{\prime 2}\mathrm{d}{r}^{\prime }.$$

(34)

The central SMBH was modelled as a point mass with values of $M_{•}=4.5_{-0.4}^{+0.9}\times {10}^9$ $M_{\odot }$, $M_{•}=6.2_{-0.5}^{+0.4}\times {10}^9$ $M_{\odot }$, and $M_{•}=1.1_{-0.1}^{+0.1}\times {10}^9$ $M_{\odot }$ for NGC 1407, NGC 4486, and NGC 5846, respectively, obtained from [64]. Therefore, the cumulative mass profile of the SMBH $M_{•}(<r)$ for each galaxy is a constant given by:

$$M_{•}(<r)=M_{•}.$$

(35)

4.2. Dynamical model

Cesare, et al. [49] modelled the RMS velocity dispersion profiles of each kinematic tracer, stars, blue GCs, and red GC, with this solution of the spherical Jeans equations [58,65,66,67]:

$$V_{\mathrm{rms},\mathrm{t}}^2\left(R\right)=\frac{2}{I_{\mathrm{t}}\left(R\right)}{\int }_R^{+\infty }K\left({\beta }_{\mathrm{t}},\frac{r}{R}\right){\nu }_{\mathrm{t}}\left(r\right)\frac{\mathrm{d}\varphi }{\mathrm{d}r}r\mathrm{d}r,$$

(36)

where (1) $\mathrm{t}=*,\mathrm{B},\mathrm{R}$ indicates the tracer, (2) $I_{\mathrm{t}}\left(R\right)$ is either the surface brightness of the stars (Equations (21) or (26)) or the 2D number density of the GCs (Equations (29) or (30)), (3) ${\nu }_{\mathrm{t}}\left(r\right)$ is either the 3D luminosity density of the stars (Equations (24) or (27)) or the 3D number density of the GCs (Equation (31)), (4) $\varphi$ is the gravitational potential and, thus, $\mathrm{d}\varphi /\mathrm{d}r$ is the gravitational field (Equation (4)), (5) R is the 2D radius (Equation (23)), (6) r is the 3D spherical radius, (7) ${\beta }_{\mathrm{t}}=1-{\sigma }_{\theta }^2/{\sigma }_r^2$ is the orbital anysotropy parameter, where ${\sigma }_{\theta }$ and ${\sigma }_r$ are the tangential and the radial velocity dispersions, and (8) K is the kernel:

$$\begin{array}{cc}\hfill K\left({\beta }_{\mathrm{t}},\frac{r}{R}\right)=& \frac{1}{2}{\left(\frac{r}{R}\right)}^{2{\beta }_{\mathrm{t}}-1}[\left(\frac{3}{2}-{\beta }_{\mathrm{t}}\right)\sqrt{\pi }\frac{\Gamma ({\beta }_{\mathrm{t}}-\frac{1}{2})}{\Gamma \left({\beta }_{\mathrm{t}}\right)}\hfill \\ \hfill & +{\beta }_{\mathrm{t}}{B}_{\frac{R^2}{r^2}}\left({\beta }_{\mathrm{t}}+\frac{1}{2},\frac{1}{2}\right)-B_{\frac{R^2}{r^2}}\left({\beta }_{\mathrm{t}}-\frac{1}{2},\frac{1}{2}\right)],\hfill \end{array}$$

(37)

given by Equation (A16) of [66], where $\Gamma \left(z\right)$ is the Euler gamma function and $B_x(a,b)$ is the incomplete beta function. The orbital anisotropy parameters ${\beta }_{*}$, ${\beta }_{\mathrm{B}}$, ${\beta }_{\mathrm{R}}$ are three of the seven free parameters of the dynamical model and they are assumed to be constant with radius [68,69].

The gravitational field $\mathrm{d}\varphi /\mathrm{d}r$ in Equation (36) is given by RG field in a spherical configuration (Equation (4)), and, thus, Equation (36) transforms into:

$$V_{\mathrm{rms},\mathrm{t}}^2\left(R\right)=\frac{2G}{I_{\mathrm{t}}\left(R\right)}{\int }_R^{+\infty }K\left({\beta }_{\mathrm{t}},\frac{r}{R}\right){\nu }_{\mathrm{t}}\left(r\right)\frac{M(<r)}{ϵ\left(\rho \right)}\frac{\mathrm{d}r}{r}.$$

(38)

The cumulative mass profile $M(<r)$ in Equation (38) is given by the sum of the cumulative mass profiles of the stars $M_{*}(<r)$ (Equation (28)), of the gas $M_{\mathrm{gas}}(<r)$ (Equation (34)) and of the SMBH $M_{•}(<r)$ (Equation (35)):

$$M(<r)=M_{*}(<r)+M_{\mathrm{gas}}(<r)+M_{•}(<r).$$

(39)

This dynamical model has seven free parameters. Four are common to the three tracers, the stellar mass-to-light ratio $\mathrm{{\rm Y}}$, and the three RG parameters, ${ϵ}_0$, Q, and ${\rho }_{\mathrm{c}}$, entering the gravitational permittivity $ϵ\left(\rho \right)$ in Equation (38), which generate the RG gravitational potential to which every kinematic tracer is subject. Three are specific to each tracer, and they are the orbital anisotropy parameters ${\beta }_{*}$, ${\beta }_{\mathrm{B}}$, and ${\beta }_{\mathrm{R}}$.

To explore the parameter space, Cesare, et al. [49] adopted a MCMC method based on Metropolis-Hastings acceptance criterion. They assumed a uniform prior on every free parameter. For $\mathrm{{\rm Y}}$, the uniform prior is based on the SPS models of Humphrey et al. [70] and Zhang et al. [61], for NGC 1407, and of Bell et al. [71] and Zibetti et al. [72], 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}\left[{\rho }_{\mathrm{c}}(\mathrm{g}/{\mathrm{cm}}^3)\right]$, respectively. At last, they explored the parameter space of the parameters ${\mathcal{B}}_{*}=-{log}_{10}(1-{\beta }_{*})$, ${\mathcal{B}}_{\mathrm{B}}=-{log}_{10}(1-{\beta }_{\mathrm{B}})$, ${\mathcal{B}}_{\mathrm{R}}=-{log}_{10}(1-{\beta }_{\mathrm{R}})$ in the uniform range of $[-1.5,1.0]$, which spans from very tangential to very radial orbits.

The model given by Equation (38) was compared with the kinematic measurements of [58], for NGC 1407, and of the ATLAS^3D survey [73], for NGC 4486 and NGC 5846, to constrain the free parameters.

4.3. Results

The RG model (38) properly interpolates the kinematic profiles of the three tracers in each E0 galaxy with sensible parameters (Figure 6).

The resulting mass-to-light ratios are consistent with the corresponding SPS models and the orbital anisotropy parameters of NGC 1407 are in agreement with the literature [58]. The RG parameters are consistent among the three galaxies within 2$\sigma$, which supports their universality (orange squares with error bars in Figure 4). The mean of the values resulting from the three single galaxies are ${ϵ}_{0,\mathrm{Mean},\mathrm{Ell}}=0.{089}_{-0.035}^{+0.038}$, $Q_{\mathrm{Mean},\mathrm{Ell}}=0.{47}_{-0.21}^{+0.29}$, and ${log}_{10}{\left[{\rho }_{\mathrm{c}}(\mathrm{g}/{\mathrm{cm}}^3)\right]}_{\mathrm{Mean},\mathrm{Ell}}=-24.{25}_{-0.20}^{+0.28}$. To better check their universality, Cesare, et al. [49] compared these values with the mean and the unique combinations of RG parameters resulting from the DMS, ${ϵ}_{0,\mathrm{Mean},\mathrm{DMS}}=0.56\pm 0.16$, $Q_{\mathrm{Mean},\mathrm{DMS}}=0.92\pm 0.71$, ${log}_{10}{\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]}_{\mathrm{Mean},\mathrm{DMS}}=-25.30\pm 1.22$ (purple squares with error bars in Figure 4), and ${ϵ}_{0,\mathrm{Unique},\mathrm{DMS}}=0.{661}_{-0.007}^{+0.007}$, $Q_{\mathrm{Unique},\mathrm{DMS}}=1.{79}_{-0.26}^{+0.14}$, ${log}_{10}{\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]}_{\mathrm{Unique},\mathrm{DMS}}=-24.{54}_{-0.07}^{+0.08}$ (green dots with blue shaded posterior distributions in Figure 4). The RG parameters derived from the elliptical galaxies are consistent with the mean combination obtained in the DMS: Q and ${log}_{10}\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]$ within 1$\sigma$ and ${ϵ}_0$ within 3$\sigma$. The parameters Q and ${log}_{10}\left[{\rho }_{\mathrm{c}}(g/c{m}^3)\right]$ from the E0 galaxies are also in agreement with the unique values from the DMS, within $3.4$$\sigma$ and $1.4$$\sigma$, respectively. Instead, the ${ϵ}_0$ parameter presents a $14.8$$\sigma$ tension.

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, $\mathrm{{\rm Y}}$ 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 said in Section 3.3, a more complete analysis could be performed by exploring the 63 parameter space of the $2\times 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 $ϵ\left(\rho \right)$, or, a fundamental problem of RG. One of the approximations introduced in the analysis of elliptical galaxies was to consider these systems as isolated, whereas these kind 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 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 [49]. Moreover, or in alternative, the functional form of the gravitational permittivity (Equation (7)) 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 [29]. If all these reasons will not be able to alleviate 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 the elliptical and the disk galaxies is obtained for the ${\rho }_{\mathrm{c}}$ parameter, which sets the density scale where the transition between the Newtonian and the RG regimes occurs, as $a_0$ in MOND.

5. 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 galaxy scale. Therefore, a class of tests at larger scale has been needed to be performed for a more complete investigation of the theory.

Some preliminary studies have already been performed by Matsakos and Diaferio [28] to study the viability of RG in modelling the hot X-ray emitting gas temperature radial profile of galaxy clusters. Matsakos and Diaferio [28] 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. These data are taken with the Chandra satellite [74]. The model that Matsakos and Diaferio [28] 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 (7) 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 ${\rho }_{\mathrm{c}}={10}^{-24}$ g cm⁻³. At last, they assumed in A1991 the presence of a stellar component of mass $5\times {10}^{10}$ M_⊙ up to a distance from the cluster centre of $R<10$ kpc. This model resulted in agreement with the Chandra data of the two galaxy clusters.

This model is certainly simplistic and it involves a too small sample of galaxy clusters. It, therefore, 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 perfomed 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 galaxy scale, to better test their universality.

Another aspect of RG on cluster scale that Matsakos and Diaferio [28] 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. 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 gravity is that it is enhanced by ${ϵ}_0^{-1}$. With these assumptions, the following expression for the BTFR would hold for galaxy groups and clusters:

$${\sigma }_V^4=GAf\frac{km}{{ϵ}_0},$$

(40)

where ${\sigma }_V$ is the velocity dispersion of the galaxies in the cluster, k is the number of galaxies, m is the mass of each galaxy, A 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 $\propto R^{-1}$. 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 ${\rho }_{\mathrm{c}}$ and to constrain ${ϵ}_0$.

Figure 8. The BTFR up to mass and velocity scales to include galaxy groups and galaxy clusters, besides galaxies. The solid line is the BTFR given by Equation (6) and the dashed line is the BTFR given by Equation (40) for $f/{ϵ}_0=1.3$. The open circles and squares represent simulated galaxies [28]. The black dots are observational data from [75,76,77]. The figure is reproduced from Figure 17 in [28].

Figure 8. The BTFR up to mass and velocity scales to include galaxy groups and galaxy clusters, besides galaxies. The solid line is the BTFR given by Equation (6) and the dashed line is the BTFR given by Equation (40) for $f/{ϵ}_0=1.3$. The open circles and squares represent simulated galaxies [28]. The black dots are observational data from [75,76,77]. The figure is reproduced from Figure 17 in [28].

6. Covariant Refracted Gravity

A covariant version of RG was recently formulated by Sanna, et al. [29]. The fact that in RG the transition between the Newtonian and the modified gravity regimes is regulated by a scalar quantity, the mass density $\rho$, rather than by a vector quantity, such as the acceleration as in MOND, made the built of a covariant extension of RG less challenging than in MOND (e.g. [78]). Some relativistic extensions of MOND, such as Tensor Vector Scalar gravity (TeVeS) [78], were formulated while presenting some problems [79,80], even if more recent results might look more promising [81,82].

Covariant refracted gravity (CRG) is formulated as a scalar tensor theory with the presence of a single scalar field, $\phi$, nonminimally coupled to the metric, which accounts for the phenomenologies both of DM on galaxy scale and of 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., [83,84,85,86,87,88,89,90,91,92,93,94,95]), and is suggested by the $a_0\sim H_0\sim \sqrt{\mathrm{\Lambda }}$ observed intriguing coincidence [29] (see Section 1).

The general action of scalar-tensor theories is:

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

(41)

where the functional forms of the self-interaction potential $\mathcal{V}\left(\phi \right)$ and the general differentiable function $\mathcal{W}\left(\phi \right)$ of the scalar field $\phi$ define a specific scalar-tensor theory. For CRG, $\mathcal{V}\left(\phi \right)$ and $\mathcal{W}\left(\phi \right)$ are:

$$\mathcal{V}\left(\phi \right)=-\Xi \phi ,$$

(42)

and

$$\mathcal{W}\left(\phi \right)=-1,$$

(43)

where $\Xi$ is a constant. Replacing Equations (42) and (43) in Equation (41), the CRG action becomes:

$${\mathcal{S}}_{\mathrm{CRG}}=\frac{1}{16\pi G}\int {\mathrm{d}}^{4}x\sqrt{g}\left[\phi R-\frac{{\nabla }^{\alpha }\phi {\nabla }_{\alpha }\phi }{\phi }-2\Xi \phi \right]+\int {\mathrm{d}}^{4}x\sqrt{g}{\mathcal{L}}_{\mathrm{m}}(g_{\mu \nu },{\psi }_{\mathrm{m}}).$$

(44)

In the WFL of the theory, the scalar field $\phi$ results to be twice the gravitational permittivity $ϵ\left(\rho \right)$. Indeed, the WFL of CRG yields the original classical formulation of RG:

$$\nabla ·(\phi \nabla \varphi )\simeq 8\pi G\rho ,$$

(45)

which is RG Poisson equation (Equation (1)) if $\phi =2ϵ\left(\rho \right)$. This result confirms that the scalar field mimics the phenomenology of DM on 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 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 ${\rho }_{\mathrm{s}}\left(r\right)$, monotonically decreasing with r, settling in a homogeneous background with constant density ${\rho }_{\mathrm{bg}}$. Sanna, et al. [29] found that, far from this source, the transition from Newtonian to RG regimes occurs when the acceleration $\mathrm{d}\varphi /\mathrm{d}r$ goes below the acceleration scale:

$$a_{\Xi }=\sqrt{2\Xi -\frac{8\pi G\rho }{\phi }},$$

(46)

where $\rho \left(r\right)={\rho }_{\mathrm{s}}\left(r\right)+{\rho }_{\mathrm{bg}}$. This recalls MOND theory, where the acceleration scale $a_0$ regulates the gravity behaviour. At large distances from the spherical source, the following limit holds:

$$2\Xi \gg 8\pi G\rho /\phi ,$$

(47)

and, therefore, Equation (46) reduces to:

$$a_{\Xi }\sim \sqrt{2\Xi }.$$

(48)

Being $\Xi \sim \mathrm{\Lambda }$, as found by an independent calculation [29] explained below, Equation (48) implies that $a_{\Xi }\sim {10}^{-10}$ m s⁻², which coincides with MOND acceleration scale $a_0$.

Sanna, et al. [29] also derived with CRG the modified Friedmann equations for a flat, expanding, homogeneous, and isotropic Universe described with the Friedmann-Lema$\widehat{i}$tre-Robertson-Walker (FLRW) metric, where the Universe content is modelled as a perfect fluid. In these modified Friedmann equations, the term

$${\Omega }_{\Xi }=\frac{\Xi }{3H^2}$$

(49)

appears, where $H\left(t\right)$ is the Hubble parameter. This term is analogous to the density parameter related to the cosmological constant $\mathrm{\Lambda }$ in $\mathrm{\Lambda }$CDM:

$${\Omega }_{\mathrm{\Lambda }}=\frac{\mathrm{\Lambda }}{3H^2},$$

(50)

which indicates that $\Xi \sim \mathrm{\Lambda }$, i.e, that $\Xi$ plays the role of $\mathrm{\Lambda }$ in $\mathrm{\Lambda }$CDM, accounting for the accelerated expansion of the Universe. At the present time, Equations (49) and (50) become:

$${\Omega }_{\Xi ,0}=\frac{\Xi }{3H_0^2},$$

(51)

and

$${\Omega }_{\mathrm{\Lambda },0}=\frac{\mathrm{\Lambda }}{3H_0^2},$$

(52)

where $H_0$ is the Hubble parameter at the present time. We have seen before that $\Xi$ sets the value of the acceleration scale $a_{\Xi }$ (Equation (46)) 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. Being $\Xi \sim \mathrm{\Lambda }$, $\Xi$ also accounts for DE phenomenology, providing a unification of the two dark sectors. Moreover, inserting $\Xi$ in Equation (46), the observed relation $a_0\sim \sqrt{\mathrm{\Lambda }}$ naturally emerges in CRG.

By rearranging the modified Friedmann equations in CRG, Sanna, et al. [29] found two solutions for the Hubble parameter $H\left(t\right)$, that refer as CRG- and CRG+ and, properly integrating $H\left(t\right)$, they found the corresponding two solutions for the scale factor $a\left(t\right)$. Properly combining CRG- and CRG+, the bound $\Xi >0$ is derived. From the derived $a\left(t\right)$, the luminosity distance $D_{\mathrm{L}}$ can be calculated and constrained from the observed Hubble diagram of SNaeIa at high redshift z of the Supernova Cosmology Project Union $2.1$ Compilation [96]. Comparing the CRG model with the data, the cosmological parameters can be derived to verify if the tensions observed in $\mathrm{\Lambda }$CDM can be either reduced or cancelled, which would provide a fundamental test for CRG. By comparing the distance modulus $\mu =m-M$, where m and M are the apparent and absolute magnitudes, respectively, derived from the luminosity distance $D_{\mathrm{L}}$ computed in CRG for the solutions CRG- and CRG+ with the $(z,\mu )$ data of the high-redshift SNaeIa from [96], Sanna, et al. [29] found that ${\Omega }_{\Xi ,0}=0.{650}_{-0.650}^{+0.005}$, for CRG-, and ${\Omega }_{\Xi ,0}=0.{650}_{-0.085}^{+0.005}$, for CRG+, assuming the values $H_0=67.7$ km s⁻¹ Mpc⁻¹ and ${\Omega }_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 ${\Omega }_{\Xi ,0}=0.65$ and $H_0=67.7$ km s⁻¹ Mpc⁻¹ in Equation (52) implies again that $\Xi \sim \mathrm{\Lambda }$, as independently determined before.

Figure 9. Hubble diagram of SNaeIa 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⁻¹ Mpc⁻¹, ${\Omega }_0=0.31$, and ${\Omega }_{\Xi ,0}=0.65$ parameters are adopted. The open circles with error bars are the data from the Supernova Cosmology Project Union 2.1 Compilation [96]. The figure is reproduced from Figure E.2 in [29].

Figure 9. Hubble diagram of SNaeIa 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⁻¹ Mpc⁻¹, ${\Omega }_0=0.31$, and ${\Omega }_{\Xi ,0}=0.65$ parameters are adopted. The open circles with error bars are the data from the Supernova Cosmology Project Union 2.1 Compilation [96]. The figure is reproduced from Figure E.2 in [29].

Sanna, et al. [29] derived the equation of state of the effective DE in CRG, $w_{\mathrm{DE}}=p_{\mathrm{DE}}/{\rho }_{\mathrm{DE}}$, where $p_{\mathrm{DE}}$ and ${\rho }_{\mathrm{DE}}$ are the pressure and the density of the effective DE, by properly rearranging the CRG Friedmann equations and by comparing them with the Friedmann equations of a general scalar-tensor theory with a non-minimal coupling between the scalar field and the metric. The $w_{\mathrm{DE}}$ parameter depends on redshift z. At the present time, $z=0$, the $w_{\mathrm{DE}}$ parameter becomes:

$$w_{\mathrm{DE}}=-\frac{6\pm 4\sqrt{3}\sqrt{1+{\Omega }_0-2{\Omega }_{\Xi 0}}}{6-3{\Omega }_0}.$$

(53)

Inserting $w_{\mathrm{DE}}=-1$, consistent with the equation of state of DE in $\mathrm{\Lambda }$CDM [1], and ${\Omega }_0\sim 0.3$, ${\Omega }_{\Xi 0}\sim 0.64$ is obtained, which is consistent with the CRG- solution.

Values of $w_{\mathrm{DE}}$ different from -1 are anyway consistent with several pieces of evidence (e.g. [97,98,99,100,101,102]). 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 [102]. For CRG, the parametrisation [103,104]:

$$w_{\mathrm{DE}}\left(z\right)=w_0+\frac{w_{a}z}{1+z}.$$

(54)

can be assumed for $w_{\mathrm{DE}}$. At the present time ($z=0$), Equation (54) only depends on $w_0$, allowing a wide range of DE models, either with $w_{\mathrm{DE}}<-1$ (phantom models) or with $w_{\mathrm{DE}}\ge -1$ [97,98,99,105,106]. The parameter $w_0$ is constrained to be approximately in the range of $w_0\in [-1.18,-0.85]$ by measurements of the baryonic acoustic oscillation (BAO), SNaeIa, and CMB [100,107]. This limit on $w_0$ translates in a limit on ${\Omega }_{\Xi 0}$, in agreement with the values of ${\Omega }_{\Xi 0}$ estimated from the SNaeIa data.

7. Discussion and conclusions

In this work, the main analyses and results of the theory of modified gravity Refracted Gravity (RG) are summarised. RG demonstrated to be able 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 galaxy scale, and to explain both the DM and the DE sectors with a single scalar field.

Several further studies could be performed to complete the tests of this gravity theory. On 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 (19) and some galaxy properties, and (2) the tension presented by the vacuum permittivity ${ϵ}_0$, which might indicate its nonuniversality.

The RAR built with RG presents some strong correlations, at more than 5$\sigma$, with some radially-dependent galaxy properties. However, also the RAR built from DMS data shows some, despite 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 DMS data are apparently at odds with the claimed uncorrelations observed in the RAR derived from SPARC data [50], 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 is observed. Di Paolo, et al. [108] 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 [50], that 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, 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 (19) of the obtained RAR have to be compared with the results of Lelli et al. [52]. Another interesting study that can 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. [108]. This study might shed light on the works of Santos-Santos, et al. [109] and Di Paolo, et al. [108], which found that the RAR might differently behave in these two categories of galaxies, besides providing a further test for RG theory.

As said at the end of Section 4.3, the tension of the vacuum permittivity ${ϵ}_0$ might have different origins, either the approximate derivation of the unique ${ϵ}_0$ in the DMS sample, or a too simplistic dynamical model for the elliptical galaxies, or an incorrect assumption for the functional form of the gravitational permittivity $ϵ\left(\rho \right)$ (Equation (7)), or, in the worst case, a fundamental issue in 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 ths new study. The E0 galaxies in the ePN.S survey [57], where extended kinematic profiles up to ∼13 $R_{\mathrm{e}}$ from the galaxy centres are 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 [110,111,112]. 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 ${\rho }_{\mathrm{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 [110,111,112] would show whether they can be reproduced by a $\propto R^{-2}$ field in the galaxy outskirts or if they necessarily require a $\propto 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 to investigate the positive correlation found by Deur [30,31] 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 [30,31] might be crucial to validate or reject this correlation. At last, all these studies might shed light on the correctness of Equation (7) for the gravitational permittivity $ϵ\left(\rho \right)$.

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 dwarf and LSB galaxies and GCs, besides of 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 galaxy scale and the latter being nearly DM-free (Section 1), which represented a challenge for MOND (e.g., [32,113,114,115,116,117,118]). 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 [119], in the LITTLE THINGS survey [120] in SPARC [50], and in Di Paolo, et al. [108], and the internal velocity dispersion of a sample of GCs located in the most external regions of the Milky Way (e.g., [32,113,114,115,116,117,118]), where the background density drops below ${\rho }_{\mathrm{c}}$.

On the scale of galaxy clusters, the tests are incomplete and a deeper investigation is required. The sample considered by Matsakos and Diaferio [28] is too small, counting two galaxy clusters, and the analysis of the temperature profiles is approximated. For their study, Matsakos and Diaferio [28] considered the data from [74], whose sample contains other 11 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 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.

Besides on 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 and HeCS, 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 to investigate how the RG parameters are affected by the environment, either relaxed or interacting, extending on larger scales the possible future study for the elliptical galaxies. Also the study of the BTFR extended to galaxy groups and clusters (Equation (40)) 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.

Also CRG leaves room for additional investigation. The emergence of MOND acceleration scale $a_0$ in the WFL of CRG should be further explored. Indeed, differently from $a_0$, $a_{\Xi }$ (Equation (46)) 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 ${\rho }_{\mathrm{s}}\ll {\rho }_{\mathrm{bg}}$. Future studies are required to test if the connection between $a_{\Xi }$ and $a_0$ subsists for a generic case besides to a source with a specific density configuration.

The parameter $w_{\mathrm{DE}}=p_{\mathrm{DE}}/{\rho }_{\mathrm{DE}}$ of the equation of state of the effective DE found by Sanna, et al. [29] in CRG appears to be dependent on redshift, differently from $\mathrm{\Lambda }$CDM where $w_{\mathrm{DE}}$ is a constant equal to -1. However, $w_{\mathrm{DE}}$ tends to -1 at the present epoch, in agreement with the observations. The $w_{\mathrm{DE}}$ parameter in CRG could be constrained from the measurements of the expansion rate of the Universe performed by the Dark Energy Survey (DES)2 [121], which have observed thousands of supernovae since August 31^th 2013. Moreover, the Euclid3 mission [122], launched on July 1^st 2023, is expected to measure the variation of the cosmic acceleration to an accuracy better than the 10% level, which would be crucial to disentangle CRG from $\mathrm{\Lambda }$CDM and other DE models [123]. 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 (54) and, consequently, the value of $\Xi$ (see Section 6).

²https://www.darkenergysurvey.org

³https://sci.esa.int/web/euclid

The Lagrangian density of CRG, derived from the CRG action (Equation (44)) 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. [124,125]), for example, by comparing the CRG predictions with the latest measurements from the Planck satellite [126]. By estimating the cosmological parameters from different low-redshift galaxy surveys, such as Kilo-Degree Survey (KIDS) 4 [127], CFHTLS 5 [128], 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 $\mathrm{\Lambda }$CDM model.

⁴https://kids.strw.leidenuniv.nl

⁵https://www.cfht.hawaii.edu/Science/CFHLS/

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., [129,130,131,132]). A possible way to accomplish this task, would be to modify the publicly available code GADGET-2 [133], 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) 6 [134], which started to observe in 2019, with Euclid, and with upcoming experiments such as the Square Kilometer Array (SKA) 7 [135].

⁶https://www.desi.lbl.gov

⁷https://www.skao.int/