###### Abstract

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.

###### keywords:

modified gravity; dark matter; galaxy dynamics; acceleration scale; scaling
relations; galaxy clusters; cosmology

3 2 6 Academic Editor: Paolo Salucci 31 January 2024 1 March 2024 1 April 2024
5 April 2024 https://doi.org/10.3390/astronomy3020006 Refracted Gravity
Solutions from Small to Large Scales Refracted Gravity Solutions from Small to
Large Scales Valentina Cesare Valentina Cesare Cesare, V. Correspondence:
<EMAIL_ADDRESS>
##  1 Introduction

In the Universe, several pieces of evidence cannot be reconciled with general
relativity (GR), if we only assume the presence of visible (baryonic) matter.
A mass discrepancy of $\sim$80–90% is observed from the largest scales (e.g.,
the cosmic microwave background (CMB) Planck Collaboration et al. (2020), the
large-scale structure Davis et al. (1985); Springel et al. (2006),
gravitational lensing effects where the Bullet Cluster example is worth
mentioning Markevitch (2006); Clowe et al. (2006); Paraficz et al. (2016),
light elements’ abundances Cyburt et al. (2016), and galaxy clusters dynamics
Zwicky (1933)) to the smallest scales (e.g., the flat trend of disk galaxies’
rotation curves Rubin and Ford (1970); Bosma (1978); Sanders (1990)).

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) Dodelson et al. (1996).
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 $\sim$70%. The most studied cosmological model is the
$\Lambda$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 $\sim$70% and $\sim$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) Kirshner (1996); Peebles and Ratra
(2003), and it can be identified with the cosmological constant $\Lambda$.

Even if $\Lambda$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 Weinberg (1989); Luo
(2014), the coincidence problem Velten et al. (2014), and the tensions between
the values of some cosmological parameters measured from probes of the late
and early Universe Planck Collaboration et al. (2020); Fleury et al. (2013);
Douspis et al. (2019). On small scales, we observe several discrepancies
between $\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 Del Popolo and Le Delliou (2017); de Martino et
al. (2020) for a review). On a galaxy scale, another issue related to the
dynamics of galaxies is present in $\Lambda$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 van Albada
and Sancisi (1986); Sackett (1997); Courteau and Rix (1998); Bissantz and
Gerhard (2002); Sellwood and Debattista (2014); McGaugh and Schombert (2015),
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 Dubinski and Carlberg (1991); Navarro et al. (1996). 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 $\sigma\sim 10$ km s-1, an
order of magnitude larger than the velocity dispersion $\sigma\sim 1$ km s-1
for systems with the same luminosity and scale radius ($\sim$100 pc) at
equilibrium Mateo (1998). 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., Torrealba et al. (2019)), which does not naturally emerge in
$\Lambda$CDM cosmological simulations and gives rise to the cusp/core problem
(e.g., Del Popolo and Le Delliou (2017); de Martino et al. (2020)). Instead,
there are systems with baryonic masses similar to those of dwarf galaxies, the
globular clusters (GCs), which are nearly DM-free Baumgardt et al. (2005,
2009); Jordi et al. (2009); Sollima and Nipoti (2010); Ibata et al. (2011a,
b); Frank et al. (2012).

It is remarkable to cite also the presence of some observed regularities on a
galaxy scale, which are hard to explain in a $\Lambda$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 $\sim$90% of the galaxy content, whereas baryonic matter only
represents $\sim$10%. Among these relations, we can mention (1) the baryonic
Tully–Fisher relation (BTFR) McGaugh et al. (2000), (2) the mass
discrepancy–acceleration relation (MDAR) McGaugh (2004), and (3) the radial
acceleration relation (RAR) McGaugh et al. (2016). The three relations see the
emergence of the same acceleration scale $a_{0}\simeq 1.2\times 10^{-10}$ m
s-2, which, if expressed in natural units, is $a_{0}\sim
H_{0}\sim\sqrt{\Lambda}$, suggesting a connection between the DM and DE
sectors. These coincidences are even less intuitive in a $\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 way is modified Newtonian dynamics (MOND)
Milgrom (1983a, b, c). 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 Matsakos and Diaferio (2016) 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, $\epsilon(\rho)$, 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 Sanna et
al. (2023) 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.

##  2 Refracted Gravity

RG is a theory of modified gravity formulated in a classical way, which can be
interpreted in analogy to electrodynamics in matter Matsakos and Diaferio
(2016): 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 Matsakos
and Diaferio (2016)):

$\nabla\cdot\left[\epsilon(\rho)\nabla\phi\right]=4\pi{G}\rho,$ (1)

where $\phi$ is the gravitational potential, and $\epsilon(\rho)$ is the
gravitational permittivity. The following asymptotic limits are adopted for
the gravitational permittivity (Equation (2.5) in Matsakos and Diaferio
(2016)):

$\epsilon(\rho)=\begin{cases}1,&\rho\gg\rho_{\rm c}\\\
\epsilon_{0},&\rho\ll\rho_{\rm c},\end{cases}$ (2)

where $0<\epsilon_{0}\leq 1$, and $\rho_{\rm 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_{\rm c}$, the Newtonian
Poisson equation is recovered (Equation (2.7) in Matsakos and Diaferio
(2016)):

$\nabla^{2}\phi_{\rm N}=4\pi{G}\rho,$ (3)

whereas, when $\rho\ll\rho_{\rm 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 Cesare et al. (2022)):

$\frac{\partial\phi}{\partial
r}=\frac{1}{\epsilon(\rho)}\frac{GM(<r)}{r^{2}}=\frac{1}{\epsilon(\rho)}\frac{\partial\phi_{\rm
N}}{\partial r},$ (4)

where $M(<r)$ is the mass of the spherical system enclosed within the
spherical radius $r$, and $\partial\phi_{\rm 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_{\rm c}$, i.e., where
$\epsilon(\rho)\rightarrow\epsilon_{0}$ (right panels of Figure 1). If we
consider a solid sphere immersed in vacuum ($\rho=0$) with radius $r_{\rm s}$
and uniform density $\rho_{\rm s}\gg\rho_{\rm c}$, the RG gravitational field
inside and outside the sphere is given by Matsakos and Diaferio (2016):

$g=\begin{cases}g_{\rm N},&r<r_{\rm s}\\\ \frac{g_{\rm
N}}{\epsilon_{0}},&r>r_{\rm s},\end{cases}$ (5)

where $g=\partial\phi/\partial{R}$ and $g_{\rm N}=\partial\phi_{\rm
N}/\partial{R}$ are the RG and Newtonian gravitational fields, respectively.
The gravitational field has the Newtonian direction everywhere but is enhanced
by $1/\epsilon_{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
Cesare et al. (2020)):

$\frac{\partial\epsilon}{\partial\rho}\nabla\rho\cdot\nabla\phi+\epsilon(\rho)\nabla^{2}\phi=4\pi
G\rho.$ (6)

In this configuration, the RG field depends both on the density $\rho$ (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
“$\frac{\partial\epsilon}{\partial\rho}\nabla\rho\cdot\nabla\phi$” 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 Deur (2014, 2020) and the different DM quantity of GCs
and dwarf galaxies, where the former are nearly spherical and DM-free (e.g.,
Sollima and Nipoti (2010)) and the latter are flatter and one of the darkest
objects in the Universe (e.g., Mateo (1998)) (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.

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 Matsakos and Diaferio (2016).

Where $\rho\ll\rho_{\rm c}$ in flat systems, e.g., in the outskirts of disk
galaxies, the radial gravitational field assumes the asymptotic trend
(Equation (B.23) in Cesare et al. (2020)):

$\frac{\partial\phi}{\partial{R}}=\sqrt{\Big{|}\frac{\partial\phi_{\rm
N}}{\partial{R}}\Big{|}a_{0}}\propto R^{-1}$ (7)

where $a_{0}$ is the MOND critical acceleration set by the normalisation of
the BTFR McGaugh (2012):

$v_{\rm f}^{4}=Ga_{0}M,$ (8)

where $v_{\rm 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 Matsakos and Diaferio (2016) derived the
relation $v_{\rm f}^{4}=GbM$ (Matsakos and Diaferio Matsakos and Diaferio
(2016) 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 McGaugh
(2012); Famaey and McGaugh (2012). 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\ll 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 ($\propto R^{-2}$)
even where $\rho\ll\rho_{\rm c}$, preserving the Gauss theorem, but the field
trend in MOND is $\propto R^{-1}$ where $a\ll a_{0}$, as in flattened systems.

In the regions of flat systems where $\rho\ll\rho_{\rm 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 Cesare et al. (2020)):

$\frac{\partial\phi}{\partial{z}}=0,$ (9)

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 ($\rho=0$) with radius $R_{\rm d}$, height
$h_{\rm d}$ from its equatorial plane (2$h_{\rm d}$ of total height), and
uniform density $\rho_{\rm d}\gg\rho_{\rm 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=\pm h_{\rm 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 $\rho$ was employed (Equation (4.1) in Matsakos and Diaferio (2016)),
Equation (5) in Cesare et al. (2020), and Equation (4) in Cesare et al.
(2022):

$\epsilon(\rho)=\epsilon_{0}+(1-\epsilon_{0})\frac{1}{2}\left\\{\tanh\left[\text{ln}\left(\frac{\rho}{\rho_{\mathrm{c}}}\right)^{Q}\right]+1\right\\},$
(10)

represented in Figure 2. Besides $\epsilon_{0}$ and $\rho_{\rm 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.

Figure 2: Gravitational permittivity $\epsilon(\rho)$ adopted in RG (Equation
(10)) for three different values of the RG free parameter $Q$ and for
$\epsilon_{0}=0.25$. The larger the value of $Q$, the steeper the transition
between the two asymptotic limits of Equation (2). The black dashed line
represents $\epsilon_{0}=0.25$. The figure is reproduced from Figure 1 in
Cesare et al. (2020). Credit: Cesare V., Diaferio A., Matsakos T., and Angus
G., A&A, 637, A70, 2020, reproduced with permission $\copyright$ ESO.

##  3 Dynamics of Disk Galaxies

Matsakos and Diaferio Matsakos and Diaferio (2016) 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. Cesare et al. (2020)
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)
Bershady et al. (2010) 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. Cesare et al. (2020) 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 Cesare (2021).

Cesare et al. Cesare et al. (2020) 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.

###  3.1 Mass Model

Cesare et al. Cesare et al. (2020) modelled the baryonic mass density profile
$\rho$, which generates the total gravitational potential $\phi$ 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 razor-thin 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 Cesare
(2021)):

$\rho_{\rm d}(R,z)=\frac{\Upsilon{I_{\rm
d}(R)}}{2h_{z}}\exp\left(-\frac{|z|}{h_{z}}\right),$ (11)

where $I_{\rm 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
$\Upsilon$ 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
Cesare et al. (2020).

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” Sancisi (2004). The disk surface
brightness in the innermost galaxy regions, where the bulge is dominant, was
modelled with an exponential profile de Vaucouleurs (1959); Freeman (1970)
(Equation (A.1) in Cesare et al. (2020)):

$I_{\rm d}(R)=I_{{\rm d},0}\exp\left(-\frac{R}{h_{R}}\right),$ (12)

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 Cesare (2021)):

$\rho_{\rm b}(r)=-\frac{\Upsilon}{\pi}\int^{+\infty}_{r}\frac{{\rm d}I_{\rm
b}(R)}{{\rm d}R}\frac{{\rm d}R}{\sqrt{R^{2}-r^{2}}},$ (13)

where $I_{\rm b}(R)$ is the surface brightness radial profile of the bulge,
and

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

is the 3D radius. The bulge contribution being subdominant, Cesare et al.
Cesare et al. (2020) assumed the mass-to-light ratio of the bulge equal to the
mass-to-light ratio of the disk, $\Upsilon$, without introducing an additional
degree of freedom. They modelled the surface brightness of the bulge with a
Sérsic spherical profile Sérsic (1963) (Equation (A.2) in Cesare et al.
(2020)):

$I_{\rm b}(R)=I_{\rm e}\exp\left\\{-7.67\left[\left(\frac{R}{R_{\rm
e}}\right)^{1/n_{\rm s}}-1\right]\right\\}.$ (15)

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
razor-thin disks (Equation (A.14) in Cesare et al. (2020)):

$\rho_{\rm atom}(R,z)=\Sigma_{\rm atom}(R)\delta(z),$ (16)

and

$\rho_{\rm mol}(R,z)=\Sigma_{\rm mol}(R)\delta(z),$ (17)

where $\Sigma_{\rm atom}(R)$ and $\Sigma_{\rm mol}(R)$ are the surface mass
densities of the atomic and molecular gas, respectively, which Cesare et al.
Cesare et al. (2020) modelled with a linear interpolation of the data as for
the stellar disk, to catch the features of the rotation curve, and $\delta(z)$
is the Dirac $\delta$ function.

The total 3D mass density profile is given by summing Equations (11), (13),
(16) and (17) (Equation (11) in Cesare (2021)):

$\rho(R,z)=\rho_{\rm d}(R,z)+\rho_{\rm b}(r)+\rho_{\rm atom}(R,z)+\rho_{\rm
mol}(R,z).$ (18)

###  3.2 Dynamical Model

From the total 3D mass density profile given by Equation (18), Cesare et al.
Cesare et al. (2020) 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 Young (1954). From the RG
potential, they derived the models for the rotation curve, $v(R)$ (Equation
(6) in Cesare et al. (2020)),

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

and the vertical velocity dispersion radial profile, $\sigma_{z}(R)$ (Equation
(13) in Cesare et al. (2020)),

$\sigma_{z}(R)=\sqrt{\frac{1}{h_{z}}\int_{0}^{+\infty}\left[\int_{z}^{+\infty}\exp\left(-\frac{|z^{\prime}|}{h_{z}}\right)\frac{\partial\phi(R,z^{\prime})}{\partial
z^{\prime}}dz^{\prime}\right]dz},$ (20)

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
$\Upsilon$, the stellar disk scale height $h_{z}$, and the three RG
parameters, $\epsilon_{0}$, $Q$, and $\rho_{\rm 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 $\Upsilon$, they adopted a
Gaussian prior based on the stellar population synthesis (SPS) model of Bell
and de Jong Bell and de Jong (2001), keeping into account that the surface
brightness was measured in the $K$-band, where the Gaussian tail for
$\Upsilon<0$ was set to zero. For $h_{z}$, they adopted a Gaussian prior
peaking at $h_{z,{\rm 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. Bershady et al. (2010) (Equation (1) in Bershady et al.
(2010) and Equation (A.12) in Cesare et al. (2020)):

$\log_{10}\left(\frac{h_{R}}{h_{z,\text{SR}}}\right)=0.367\log_{10}\left(\frac{h_{R}}{\text{kpc}}\right)+0.708\pm
0.095\;,$ (21)

where the term $\pm 0.095$ indicates the $\sim$25% intrinsic scatter.
Estimating $h_{z,{\rm 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 $\epsilon_{0}$,
$[0.01,2.00]$ for $Q$, and $[-27.00,-23.00]$ for $\log_{10}[\rho_{\rm c}({\rm
g}/{\rm cm}^{3})]$.

###  3.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,{\rm SR}}$ (Equation
(21)). In all panels, the dashed magenta vertical lines identify the bulge
effective radius ($R_{\rm 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. Cesare et al. (2020) found mass-
to-light ratios $\Upsilon$ generally consistent with the SPS model of Bell and
de Jong Bell and de Jong (2001). In the only-rotation curve analysis, the
estimated $h_{z}$’s were generally consistent with the $h_{z,{\rm SR}}$
derived from Equation (21). Instead, in the combined analysis, the estimated
$h_{z}$’s were systematically smaller than the correspondent $h_{z,{\rm SR}}$.
Angus et al. Angus et al. (2015) found a similar result by modelling the DMS
galaxies with QUMOND Milgrom (2010), 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. Cesare et al. (2020) and by Angus et al.
Angus et al. (2015) 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 Aniyan et al. (2016). The young stellar population
distributes in a disk thinner and with a smaller vertical velocity dispersion
$\sigma_{z}$ than the old stellar population Aniyan et al. (2016). Instead,
the disk scale heights of edge-on galaxies, from which Equation (21) is
derived to determine $h_{z,{\rm SR}}$, are directly measured from near-
infrared photometry, whose signal is dominated by the old stellar population
Aniyan et al. (2016), 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. Cesare et al. (2020) computed this factor, which was
equal to $1.63$, in agreement with the work of Aniyan et al. Aniyan et al.
(2016). Cesare et al. Cesare et al. (2020) 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,{\rm SR}}$.

Figure 3: 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,{\rm SR}}$ (Equation (21)), and the dashed magenta
and green vertical lines show the bulge effective radius ($R_{\rm 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 Cesare et al. (2020). Credit: Cesare V., Diaferio
A., Matsakos T., and Angus G., A&A, 637, A70, 2020, reproduced with permission
$\copyright$ 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. Cesare et al. (2020) 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 $\epsilon_{0,{\rm
Mean,DMS}}=0.56\pm 0.16$, $Q_{\rm Mean,DMS}=0.92\pm 0.71$,
$\log_{10}[\rho_{\rm c}\text{(g/cm}^{3})]_{\rm Mean,DMS}=-25.30\pm 1.22$
(purple squares with error bars in Figure 4), whereas the unique combination
of free parameters was $\epsilon_{0,{\rm
Unique,DMS}}=0.661^{+0.007}_{-0.007}$, $Q_{\rm
Unique,DMS}=1.79^{+0.14}_{-0.26}$, $\log_{10}[\rho_{\rm
c}\text{(g/cm}^{3})]_{\rm 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$\sigma$, 2$\sigma$, and 3$\sigma$ 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 Cesare et al. (2020); Aldinucci et al. (2021).

Figure 4: Free parameters of the RG gravitational permittivity
$\epsilon(\rho)$ 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$\sigma$, 2$\sigma$, and 3$\sigma$ 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 Cesare et al.
(2022). Credit: Cesare V., Diaferio A., and Matsakos T., A&A, 657, A133, 2022,
reproduced with permission $\copyright$ ESO.

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
$\epsilon_{0}$, $Q$, and $\log_{10}[\rho_{\rm c}\text{(g/cm}^{3})]$ agreed
with the correspondent RG parameters estimated from the single galaxies within
$0.63$$\sigma$, $1.19$$\sigma$, and $0.62$$\sigma$, 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
$\Upsilon$ 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. Cesare et al. (2020) also modelled the RAR of the DMS galaxies.
McGaugh et al. McGaugh et al. (2016) fitted the RAR of 153 edge-on disk
galaxies in the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample
Lelli et al. (2016) with this relation (Equation (4) in McGaugh et al. (2016)
and Equation (20) in Cesare et al. (2020)):

$g_{\text{obs}}(R)=\frac{g_{\text{bar}}(R)}{1-\exp\left(-\sqrt{\frac{g_{\text{bar}}(R)}{g_{\dagger}}}\right)},$
(22)

where $g_{\rm obs}(R)$ is the observed centripetal acceleration derived from
the measured rotation curve $v_{\rm obs}(R)$ (Equation (19) in Cesare et al.
(2020)),

$g_{\rm obs}(R)=\frac{v^{2}_{\rm obs}(R)}{R},$ (23)

$g_{\rm bar}(R)$ is the Newtonian acceleration due to baryonic matter alone,
and $g_{\dagger}$ is the only free parameter of the model. This relation
properly interpolated the SPARC data with $g_{\dagger}=[1.20\pm 0.02{\rm\text{
}(random)}\pm 0.24{\rm\text{ }(systematic)}]\times 10^{-10}$ m s-2, 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 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. Li et al. (2018)
from SPARC galaxies) and with residuals strongly correlated, largely at more
than 5$\sigma$, 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.
Lelli et al. (2017) 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.

Figure 5: 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. McGaugh et al.
(2016) from SPARC galaxies and the black dashed line is the $g_{\rm
obs}=g_{\rm bar}$ relation, plotted as a reference. The figure is readapted
from Figure 11 in Cesare et al. (2020). Credit: Cesare V., Diaferio A.,
Matsakos T., and Angus G., A&A, 637, A70, 2020, reproduced with permission
$\copyright$ ESO.

##  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
objects, Cesare et al. Cesare et al. (2022) 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 Pota et al. (2013); Brodie
et al. (2014); Forbes et al. (2017), 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 Cesare (2021).

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 $\varepsilon=1-q$
of $0.05$, $0.14$, and $0.08$, they can be approximated as spherical systems.
Cesare et al. Cesare et al. (2022) chose these galaxies since the kinematic
information of SLUGGS galaxies was measured up to their most external regions
($\gtrsim$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 Mathews and Brighenti (2003), planetary nebulae Pulsoni et al. (2018), and
GCs Pota et al. (2013), which settle in the galaxies outskirts. Specifically,
the choice of Cesare et al. (2022) 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 Pota et al. (2013). 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.

###  4.1 Mass Model

Cesare et al. Cesare et al. (2022) modelled, at the same time and for each E0
galaxy, the root-mean-square (RMS) velocity dispersions of the stars,
concentrated within about one effective radius ($R_{\rm 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 $\sim$5, $\sim$48, and $\sim$13 times more extended than the
kinematic profiles of the stars in NGC 1407, NGC 4486, and NGC 5846. Cesare et
al. Cesare et al. (2022) 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. Cesare et al. (2022) 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 $\sim$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 Pota et al. (2015), was modelled with a Sérsic profile Sérsic
(1963) of this kind, with parameters from Pota et al. (2013) (Equation (7) in
Cesare et al. (2022)):

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

where (Equation (8) in Cesare et al. (2022))

$b_{n_{\rm s}}\thickapprox 2n_{\rm s}-\frac{1}{3}+\frac{4}{405n_{\rm
s}}+\frac{46}{25515n_{\rm s}^{2}}$ (25)

is such that half of the total luminosity of the stars is enclosed within 1
$R_{\rm e}$. The circularised radius $R$ is given for every galaxy by
(Equation (6) in Cesare et al. (2022)):

$R=\sqrt{qx^{\prime 2}+\frac{y^{\prime 2}}{q}},$ (26)

where $x^{\prime}$ and $y^{\prime}$ 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
Cesare et al. (2022)):

$\nu_{*}(r)=-\frac{1}{\pi}\int^{+\infty}_{r}\frac{{\rm d}I_{*}(R)}{{\rm
d}R}\frac{{\rm d}R}{\sqrt{R^{2}-r^{2}}},$ (27)

and integrating Equation (27) in spherical coordinates yields its cumulative
stellar luminosity profile (Equation (14) in Cesare et al. (2022)):

$L_{*}(<r)=4\pi\int^{r}_{0}\nu_{*}(r^{\prime})r^{\prime 2}{\rm d}r^{\prime}.$
(28)

For the surface brightness of the stars of NGC 4486 and NGC 5846, measured in
the $r$-band Scott et al. (2013), Cesare et al. Cesare et al. (2022) adopted a
multi-Gaussian expansion (MGE) model Emsellem et al. (1994) in a spherical
configuration, with the parameters of Scott et al. (2013) (Equation (10) in
Cesare et al. (2022)):

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

The 3D luminosity density distribution of the stars of NGC 4486 and NGC 5846
was given by (Equation (13) in Cesare et al. (2022))

$\nu_{*}(r)=\sum_{k=1}^{N}\frac{L_{k}}{(\sqrt{2\pi}\sigma_{k})^{3}}\exp\left(-\frac{r^{2}}{2\sigma^{2}_{k}}\right)$
(30)

and their cumulative luminosity stellar profile $L_{*}(<r)$ was obtained by
Equation (28), where $\nu(r^{\prime})$ is Equation (30).

The cumulative mass profile of the stars in each galaxy, $M_{*}(<r)$, was
obtained by (Equation (28) in Cesare et al. (2022)):

$M_{*}(<r)=\Upsilon{L_{*}(<r)},$ (31)

where $L_{*}(<r)$ is Equation (28), and $\Upsilon$ 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 Sérsic (1963) of this kind (Equation
(15) in Cesare et al. (2022)):

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

where Cesare et al. Cesare et al. (2022) adopted the parameters of Pota et al.
(2015) for NGC 1407 and estimated the parameters from the data of Pota et al.
(2013) with an MCMC for NGC 5846. For NGC 4486, they adopted a different
parametrisation (Equation (16) in Cesare et al. (2022)):

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

with the parameters of Strader et al. (2011). The 3D number density profiles
were obtained by integrating Equations (32) and (33) with the Abel integral
(Equation (17) in Cesare et al. (2022)):

$\nu_{\rm GC}(r)=-\frac{1}{\pi}\int^{+\infty}_{r}\frac{{\rm d}N_{\rm
GC}(R)}{{\rm d}R}\frac{{\rm d}R}{\sqrt{R^{2}-r^{2}}}.$ (34)

Cesare et al. Cesare et al. (2022) modelled the mass density profile of the
gas of NGC 1407 and NGC 5846 using the two-$\beta$ functional form (Equations
(18) and (19) in Cesare et al. (2022)):

$\rho_{\rm gas}(r)=\mu{m_{\rm H}}\sqrt{n_{\rm
g,1}^{2}\left[1+\left(\frac{r}{R_{\rm
c,1}}\right)^{2}\right]^{-3\beta_{1}}+n_{\rm
g,2}^{2}\left[1+\left(\frac{r}{R_{\rm c,2}}\right)^{2}\right]^{-3\beta_{2}}},$
(35)

with the parameters of Zhang et al. (2007) for NGC 1407, and estimating the
parameters from the data of Paggi et al. (2017) with an MCMC for NGC 5846. In
Equation (35), $\mu=0.6$ and $m_{\rm 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 Fabricant et al. (1980)
(Equation (20) in Cesare et al. (2022)):

$\rho_{\rm
gas}(r)=\frac{\rho_{0}}{(1+b^{\prime}r^{2}+c^{\prime}r^{4}+d^{\prime}r^{6})^{n^{\prime}}}.$
(36)

The cumulative mass density profile of the gas, $M_{\rm gas}(<r)$, was
obtained by integrating in spherical coordinates $\rho_{\rm gas}(r)$ from
Equation (35) for NGC 1407 and NGC 5846, and Equation (36) for NGC 4486
(Equation (14) in Cesare et al. (2022)):

$M_{\rm gas}(<r)=4\pi\int^{r}_{0}\rho_{\rm gas}(r^{\prime})r^{\prime 2}{\rm
d}r^{\prime}.$ (37)

The central SMBH was modelled as a point mass with values of
$M_{\bullet}=4.5^{+0.9}_{-0.4}\times 10^{9}$ $M_{\odot}$,
$M_{\bullet}=6.2^{+0.4}_{-0.5}\times 10^{9}$ $M_{\odot}$, and
$M_{\bullet}=1.1^{+0.1}_{-0.1}\times 10^{9}$ $M_{\odot}$ for NGC 1407, NGC
4486, and NGC 5846, obtained from Rusli et al. (2013).

###  4.2 Dynamical Model

Cesare et al. Cesare et al. (2022) modelled the RMS velocity dispersion
profiles of each kinematic tracer, stars, blue GCs, and red GC, with this
solution of the spherical Jeans equations Jeans (1915); Mamon and Łokas
(2005); Cappellari (2008); Pota et al. (2015) (Equation (23) in Cesare et al.
(2022)):

$V_{\rm rms,t}^{2}(R)=\frac{2}{I_{\rm
t}(R)}\int_{R}^{+\infty}K\left(\beta_{\rm t},\frac{r}{R}\right)\nu_{\rm
t}(r)\frac{{\rm d}\phi}{{\rm d}r}r\text{ }{\rm d}r,$ (38)

where (1) ${\rm t}={*,{\rm B},{\rm R}}$ indicates the tracer, (2) $I_{\rm
t}(R)$ is either the surface brightness of the stars (Equations (24) or (29))
or the 2D number density of the GCs (Equations (32) or (33)), (3) $\nu_{\rm
t}(r)$ is either the 3D luminosity density of the stars (Equations (27) or
(30)) or the 3D number density of the GCs (Equation (34)), (4) $\phi$ is the
gravitational potential and, thus, ${\rm d}\phi/{\rm d}r$ is the gravitational
field (Equation (4)), (5) $R$ is the 2D radius (Equation (26)), (6) $r$ is the
3D spherical radius, (7) $\beta_{\rm t}=1-\sigma_{\theta}^{2}/\sigma_{r}^{2}$
is the orbital anisotropy parameter, where $\sigma_{\theta}$ and $\sigma_{r}$
are the tangential and radial velocity dispersions, and (8) $K$ is the kernel
(Equation (A16) in Mamon and Łokas (2005) and Equation (24) in Cesare et al.
(2022)):

$\begin{split}K\left(\beta_{\rm
t},\frac{r}{R}\right)=&\frac{1}{2}\left(\frac{r}{R}\right)^{2\beta_{\rm
t}-1}\bigg{[}\left(\frac{3}{2}-\beta_{\rm
t}\right)\sqrt{\pi}\frac{\Gamma(\beta_{\rm t}-\frac{1}{2})}{\Gamma(\beta_{\rm
t})}\\\ &+\beta_{\rm t}B_{\frac{R^{2}}{r^{2}}}\left(\beta_{\rm
t}+\frac{1}{2},\frac{1}{2}\right)-B_{\frac{R^{2}}{r^{2}}}\left(\beta_{\rm
t}-\frac{1}{2},\frac{1}{2}\right)\bigg{]},\\\ \end{split}$ (39)

where $\Gamma(z)$ is the Euler gamma function, and $B_{x}(a,b)$ is the
incomplete beta function. The orbital anisotropy parameters $\beta_{*}$,
$\beta_{\rm B}$, and $\beta_{\rm R}$ are three of the seven free parameters of
the dynamical model, and they are assumed to be constant with the radius
Thomas et al. (2014); Rantala et al. (2019).

The gravitational field ${\rm d}\phi/{\rm d}r$ in Equation (38) is given by
the RG field in a spherical configuration (Equation (4)), and thus, Equation
(38) transforms into (Equation (25) in Cesare et al. (2022)):

$V_{\rm rms,t}^{2}(R)=\frac{2G}{I_{\rm
t}(R)}\int_{R}^{+\infty}K\left(\beta_{\rm t},\frac{r}{R}\right)\nu_{\rm
t}(r)\frac{M(<r)}{\epsilon(\rho)}\frac{{\rm d}r}{r}.$ (40)

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_{\rm gas}(<r)$ (Equation (37)), and the SMBH (end of Section 4.1):

$M(<r)=M_{*}(<r)+M_{\rm gas}(<r)+M_{\bullet},$ (41)

taken from Equation (26) in Cesare et al. (2022).

This dynamical model has seven free parameters. Four parameters are common to
the three tracers, the stellar mass-to-light ratio $\Upsilon$, and the three
RG parameters, $\epsilon_{0}$, $Q$, and $\rho_{\rm c}$, entering the
gravitational permittivity $\epsilon(\rho)$ 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 $\beta_{*}$, $\beta_{\rm B}$, and $\beta_{\rm R}$.

To explore the parameter space, Cesare et al. Cesare et al. (2022) adopted a
MCMC method based on a Metropolis–Hastings acceptance criterion. They assumed
a uniform prior on every free parameter. For $\Upsilon$, the uniform prior was
based on the SPS models of Humphrey et al. Humphrey et al. (2006) and Zhang et
al. Zhang et al. (2007) for NGC 1407, and of Bell et al. Bell et al. (2003)
and Zibetti et al. Zibetti et al. (2009) 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 $\epsilon_{0}$, $Q$, and
$\log_{10}[\rho_{\rm c}({\rm g}/{\rm cm}^{3})]$, respectively. At last, they
explored the parameter space of the parameters
$\mathcal{B}_{*}=-\log_{10}(1-\beta_{*})$, $\mathcal{B}_{\rm
B}=-\log_{10}(1-\beta_{\rm B})$, $\mathcal{B}_{\rm R}=-\log_{10}(1-\beta_{\rm
R})$ 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 Pota et al. (2015) for NGC 1407, and of the ATLAS3D survey Cappellari et
al. (2011) for NGC 4486 and NGC 5846, to constrain the free parameters.

###  4.3 Results

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

Figure 6: Root-mean-square velocity dispersion profiles computed with RG with
the parameters resulting from the MCMC analysis of Cesare et al. (2022),
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 $\Upsilon$ and $\beta$ 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 Cesare et al.
(2022). Credit: Cesare V., Diaferio A., and Matsakos T., A&A, 657, A133, 2022,
reproduced with permission $\copyright$ ESO.

The resulting mass-to-light ratios were consistent with the corresponding SPS
models, and the orbital anisotropy parameters of NGC 1407 were in agreement
with the literature Pota et al. (2015). The RG parameters were consistent
among the three galaxies within 2$\sigma$, which supported their universality
(orange squares with error bars in Figure 4). The mean of the values resulting
from the three single galaxies were $\epsilon_{0,{\rm
Mean,Ell}}=0.089^{+0.038}_{-0.035}$, $Q_{\rm Mean,Ell}=0.47^{+0.29}_{-0.21}$,
and $\log_{10}[\rho_{\rm c}({\rm g}/{\rm cm}^{3})]_{\rm
Mean,Ell}=-24.25^{+0.28}_{-0.20}$. To better check their universality, Cesare
et al. Cesare et al. (2022) compared these values with the mean and the unique
combinations of RG parameters resulting from the DMS galaxies,
$\epsilon_{0,{\rm Mean,DMS}}=0.56\pm 0.16$, $Q_{\rm Mean,DMS}=0.92\pm 0.71$,
$\log_{10}[\rho_{\rm c}({\rm g}/{\rm cm}^{3}))]_{\rm Mean,DMS}=-25.30\pm 1.22$
(purple squares with error bars in Figure 4), and $\epsilon_{0,{\rm
Unique,DMS}}=0.661^{+0.007}_{-0.007}$, $Q_{\rm
Unique,DMS}=1.79^{+0.14}_{-0.26}$, $\log_{10}[\rho_{\rm c}({\rm g}/{\rm
cm}^{3}))]_{\rm Unique,DMS}=-24.54^{+0.08}_{-0.07}$ (green dots with blue
shaded posterior distributions in Figure 4). The RG parameters derived from
the elliptical galaxies were consistent with the mean combination obtained in
the DMS: $Q$ and $\log_{10}[\rho_{\rm c}({\rm g}/{\rm cm}^{3}))]$ within
1$\sigma$ and $\epsilon_{0}$ within 3$\sigma$. The parameters $Q$ and
$\log_{10}[\rho_{\rm c}({\rm g}/{\rm cm}^{3}))]$ from the E0 galaxies were
also in agreement with the unique values from the DMS, within $3.4$$\sigma$
and $1.4$$\sigma$, respectively. Instead, the $\epsilon_{0}$ parameter
presented a $14.8$$\sigma$ tension.

The reasons for the $\epsilon_{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, $\Upsilon$ 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 $\epsilon_{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\times 30$ mass parameters and the
three RG parameters in the DMS.

If, instead, the posterior distribution of $\epsilon_{0}$ is not particularly
underestimated, other reasons could explain the $\epsilon_{0}$ discrepancy,
such as a too approximated model for the elliptical galaxies, a wrong
functional form for the gravitational permittivity $\epsilon(\rho)$, 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 $\epsilon_{0}$s. Another
approximation adopted in modelling the elliptical galaxies was to neglect
their net rotation, although these three galaxies are slow rotators Cesare et
al. (2022). 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 Sanna et al. (2023). If all these reasons
are not able to alleviate or cancel the $\epsilon_{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 $\rho_{\rm 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 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
Matsakos and Diaferio (2016) to study the viability of RG in modelling the hot
X-ray-emitting gas temperature radial profile of galaxy clusters. Matsakos and
Diaferio Matsakos and Diaferio (2016) 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
Vikhlinin et al. (2006). These data are taken with the Chandra satellite
Vikhlinin et al. (2006). The model that Matsakos and Diaferio Matsakos and
Diaferio (2016) 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 $\epsilon_{0}=0.045$ and
$0.065$ for A1991 and A1795, respectively, $Q=2$, and $\rho_{\rm c}=10^{-24}$
g cm-3. 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 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.

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 Vikhlinin et al. (2006).
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 Matsakos and
Diaferio (2016).

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 Matsakos
and Diaferio (2016) 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 $\epsilon_{0}^{-1}$. 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 Matsakos and Diaferio (2016)):

$v_{\rm f}^{4}=Gbf\frac{kM}{\epsilon_{0}},$ (42)

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 $\propto R^{-1}$. In Equation (42), $v_{\rm 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_{\rm 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/\epsilon_{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 McGaugh (2005), having
a mass range of $M\sim[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 Sanders (2003); McGaugh (2015), having a mass range of
$M\sim[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 $\rho_{\rm c}$ and to constrain $\epsilon_{0}$.

Figure 8: The BTFR up to mass and velocity scales to include galaxy groups and
clusters, besides galax-ies. 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/\epsilon_{0}=1.3$. The open circles and squares represent simulated
galaxies Matsakos and Diaferio (2016). The black dots that concentrate in the
bottom-left part of the plot are observational data from McGaugh (2005) and
refer to galaxies. The black dots that concentrate in the top-right part of
the plot are observational data from Sanders (2003); McGaugh (2015) 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
in Matsakos and Diaferio (2016).

##  6 Covariant Refracted Gravity

A covariant version of RG was recently formulated by Sanna et al. Sanna et al.
(2023). 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 building of a covariant extension of RG less challenging than in MOND
(e.g., Bekenstein (2004)). Some relativistic extensions of MOND, such as
tensor vector scalar gravity (TeVeS) Bekenstein (2004), were formulated while
presenting some problems Skordis (2009); Bekenstein (2011), even if more
recent results might look more promising Hernandez et al. (2019); Skordis and
Złośnik (2021).

Covariant refracted gravity (CRG) is formulated as a scalar tensor theory with
the presence of a single scalar field, $\varphi$, 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., Kunz et al. (2009); Kunz (2009); Makler
et al. (2003); Bento et al. (2003); Scherrer (2004); Bruni et al. (2013);
Cadoni et al. (2018); Sotiriou and Faraoni (2010); Sebastiani et al. (2016);
Berezhiani et al. (2017); Ferreira et al. (2019); Arbey and Coupechoux (2021);
Campigotto et al. (2019)) and is suggested by the $a_{0}\sim
H_{0}\sim\sqrt{\Lambda}$ observed intriguing coincidence Sanna et al. (2023)
(see Section 1).

The general action of scalar–tensor theories is (Equation (3) in Sanna et al.
(2023)):

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

where the functional forms of the self-interaction potential
$\mathcal{V}(\varphi)$ and the general differentiable function
$\mathcal{W}(\varphi)$ of the scalar field $\varphi$ define a specific
scalar–tensor theory. For CRG, $\mathcal{V}(\varphi)$ and
$\mathcal{W}(\varphi)$ are (Equation (7) in Sanna et al. (2023)):

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

and (Equation (6) in Sanna et al. (2023))

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

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

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

In the WFL of the theory, the scalar field $\varphi$ becomes twice the
gravitational permittivity $\epsilon(\rho)$. Indeed, the WFL of CRG yields the
original classical formulation of RG (Equation (24) in Sanna et al. (2023)):

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

which is the RG Poisson equation (Equation (1)) if $\varphi=2\epsilon(\rho)$.
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 $\rho_{\rm s}(r)$, monotonically decreasing with $r$, settling in a
homogeneous background with constant density $\rho_{\rm bg}$. Sanna et al.
Sanna et al. (2023) found that, far from this source, the transition from
Newtonian to RG regimes occurs when the acceleration ${\rm d}\phi/{\rm d}r$
goes below the acceleration scale (Equation (28) in Sanna et al. (2023)):

$a_{\Xi}=\sqrt{2\Xi-\frac{8\pi G\rho}{\varphi}},$ (48)

where $\rho(r)=\rho_{\rm s}(r)+\rho_{\rm 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 Sanna et al. (2023)):

$2\Xi\gg 8\pi G\rho/\varphi,$ (49)

and therefore, Equation (48) reduces to (Section 3.4.1 in Sanna et al.
(2023)):

$a_{\Xi}\sim\sqrt{2\Xi}.$ (50)

Since $\Xi\sim\Lambda$, as found by independent calculations Sanna et al.
(2023) explained below, Equation (50) implies that $a_{\Xi}\sim 10^{-10}$ m
s-2, which coincides with the MOND acceleration scale $a_{0}$.

Sanna et al. Sanna et al. (2023) also derived with CRG the modified Friedmann
equations for a flat, expanding, homogeneous, and isotropic Universe described
with the Friedmann–Lema$\hat{i}$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 Sanna et al. (2023))

$\Omega_{\Xi}=\frac{\Xi}{3H^{2}}$ (51)

appears, where $H(t)$ is the Hubble parameter. This term is analogous to the
density parameter related to the cosmological constant $\Lambda$ in
$\Lambda$CDM (Section 4.2 in Sanna et al. (2023)):

$\Omega_{\Lambda}=\frac{\Lambda}{3H^{2}},$ (52)

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

$\Omega_{\Xi,0}=\frac{\Xi}{3H_{0}^{2}},$ (53)

and

$\Omega_{\Lambda,0}=\frac{\Lambda}{3H_{0}^{2}},$ (54)

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 (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 $\Xi\sim\Lambda$, $\Xi$ also accounts for DE
phenomenology, providing a unification of the two dark sectors. Moreover,
inserting $\Xi$ in Equation (48), the observed relation
$a_{0}\sim\sqrt{\Lambda}$ naturally emerges in CRG.

By rearranging the modified Friedmann equations in CRG, Sanna et al. Sanna et
al. (2023) 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 $\Xi>0$ was derived. From the derived $a(t)$, the
luminosity distance $D_{\rm 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 Suzuki et al. (2012) (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
$\Lambda$CDM could 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_{\rm L}$ computed in CRG for the solutions
CRG$-$ and CRG$+$ (dashed and solid lines in Figure 9) with the $(z,\mu)$ data
of the high-redshift SNe Ia from Suzuki et al. (2012), Sanna et al. Sanna et
al. (2023) found that $\Omega_{\Xi,0}=0.650^{+0.005}_{-0.650}$ for CRG$-$, and
$\Omega_{\Xi,0}=0.650^{+0.005}_{-0.085}$ for CRG$+$, assuming the values
$H_{0}=67.7$ km s-1 Mpc-1 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 Planck
Collaboration et al. (2020). Replacing the values of $\Omega_{\Xi,0}=0.65$ and
$H_{0}=67.7$ km s-1 Mpc-1 in Equation (54) implies again that
$\Xi\sim\Lambda$, as independently determined before.

Sanna et al. Sanna et al. (2023) derived the equation of state of the
effective DE in CRG, $w_{\rm DE}=p_{\rm DE}/\rho_{\rm DE}$, where $p_{\rm DE}$
and $\rho_{\rm 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_{\rm DE}$ parameter
depends on redshift $z$. At the present time, $z=0$, the $w_{\rm DE}$
parameter becomes (Equation (73) in Sanna et al. (2023)):

$w_{\rm DE}=-\frac{6\pm 4\sqrt{3}\sqrt{1+\Omega_{0}-2\Omega_{\Xi
0}}}{6-3\Omega_{0}}.$ (55)

Inserting $w_{\rm DE}=-1$, consistent with the equation of state of DE in
$\Lambda$CDM Planck Collaboration et al. (2020), and $\Omega_{0}\sim 0.3$,
$\Omega_{\Xi 0}\sim 0.64$ is obtained, which is consistent with the
CRG$-$solution.

Values of $w_{\rm DE}$ different from $-$1 are consistent with several pieces
of evidence anyway (e.g., Amendola and Tsujikawa (2015); Copeland et al.
(2006); Frusciante and Perenon (2020); Wen et al. (2018); Capozziello et al.
(2006); Gerardi et al. (2019)). 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 Gerardi et al. (2019). For
CRG, the parametrisation Chevallier and Polarski (2001); Linder (2003)
(Section 4.3 in Sanna et al. (2023)):

$w_{\rm DE}(z)=w_{0}+\frac{w_{a}z}{1+z}.$ (56)

can be assumed for $w_{\rm DE}$. At the present time ($z=0$), Equation (56)
only depends on $w_{0}$, allowing a wide range of DE models, either with
$w_{\rm DE}<-1$ (phantom models) or with $w_{\rm DE}\geq-1$ Bean and
Melchiorri (2002); Solà and Štefančić (2005); Copeland et al. (2006); Amendola
and Tsujikawa (2015); Frusciante and Perenon (2020). 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), SNe Ia, and CMB Hazra
et al. (2015); Wen et al. (2018). This limit on $w_{0}$ translates into a
limit on $\Omega_{\Xi 0}$, in agreement with the values of $\Omega_{\Xi 0}$
estimated from the SNe Ia data.

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, $\Omega_{0}=0.31$, and
$\Omega_{\Xi,0}=0.65$ parameters were adopted. The open circles with error
bars are the data from the Supernova Cosmology Project Union $2.1$ Compilation
Suzuki et al. (2012). The figure is reproduced from Figure E.2 in Sanna et al.
(2023).

##  7 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 $\epsilon_{0}$, which might indicate its nonuniversality.

Matsakos and Diaferio Matsakos and Diaferio (2016) 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$\sigma$, 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 Lelli et al.
(2016), 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. Di Paolo et al. (2019) 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 Lelli et al. (2016), 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. Lelli et al. (2017). 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. Di Paolo et al. (2019). This study might shed
light on the works of Santos-Santos et al. Santos-Santos et al. (2020) and Di
Paolo et al. Di Paolo et al. (2019), 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
$\epsilon_{0}$ might have different origins, either the approximate derivation
of the unique $\epsilon_{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 $\epsilon(\rho)$ (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 $\epsilon_{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
Pulsoni et al. (2018), where extended kinematic profiles up to $\sim$13
$R_{\rm 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 Jiménez et al. (2013); Durazo et al. (2017, 2018).
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_{\rm 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 Jiménez et al. (2013); Durazo et al. (2017, 2018) 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 one to investigate the positive correlation found by
Deur Deur (2014, 2020) 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 $\sim$10 times more extended than the data used by Deur
(2014, 2020) 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 $\epsilon(\rho)$.

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., Baumgardt et al. (2005); Jordi et al.
(2009); Baumgardt et al. (2009); Sollima and Nipoti (2010); Ibata et al.
(2011a, b); Frank et al. (2012)). 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 Salucci et al. (2012), in the LITTLE THINGS survey Oh et al.
(2015), in SPARC Lelli et al. (2016), and in Di Paolo et al. Di Paolo et al.
(2019), and the internal velocity dispersion of a sample of GCs located in the
most external regions of the Milky Way (e.g., Baumgardt et al. (2005); Jordi
et al. (2009); Baumgardt et al. (2009); Sollima and Nipoti (2010); Ibata et
al. (2011a, b); Frank et al. (2012)), where the background density drops below
$\rho_{\rm 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 $\Lambda$CDM Mateo (1998) (See
Section 1). Moreover, Matsakos and Diaferio Matsakos and Diaferio (2016)
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 $\Lambda$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 Walker et al. (2007); Salomon et al. (2015).

On the scale of galaxy clusters, the tests are incomplete, and a deeper
investigation is required. The sample considered by Matsakos and Diaferio
Matsakos and Diaferio (2016) is too small, only counting two galaxy clusters,
and the analysis of the temperature profiles is approximated. For their study,
Matsakos and Diaferio Matsakos and Diaferio (2016) considered the data from
Vikhlinin et al. (2006), 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 Rines and Diaferio (2006) and HeCS Rines et al. (2013),
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 Gianfagna et
al. (2021). 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_{\Xi}$ (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 $\rho_{\rm s}\ll\rho_{\rm bg}$.
Future studies are required to test if the connection between $a_{\Xi}$ and
$a_{0}$ subsists for a generic case besides a source with a specific density
configuration.

The parameter $w_{\rm DE}=p_{\rm DE}/\rho_{\rm DE}$ of the equation of state
of the effective DE found by Sanna et al. Sanna et al. (2023) in CRG appears
to be dependent on the redshift, differently from $\Lambda$CDM where $w_{\rm
DE}$ is a constant equal to $-$1\. However, $w_{\rm DE}$ tends to $-$1 at the
present epoch, in agreement with the observations. The $w_{\rm 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) The Dark Energy
Survey Collaboration (2005), 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) Euclid Collaboration et al. (2022), 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 $\Lambda$CDM and other DE models Amendola et al. (2018). 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 $\Xi$ (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.,
Huterer et al. (2015); Noller and Nicola (2019)), for example, by comparing
the CRG predictions with the latest measurements from the Planck satellite
Planck Collaboration et al. (2020). 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) de Jong et
al. (2013), Canada-France-Hawaii Telescope Legacy Survey (CFHTLS)
(https://www.cfht.hawaii.edu/Science/CFHLS/, accessed on 4 April 2024) Gwyn
(2012), 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
$\Lambda$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., di Porto and Amendola (2008); Bueno Sánchez
and Perivolaropoulos (2011); Pace et al. (2014); Kofinas and Lima (2017)). A
possible way to accomplish this task, would be to modify the publicly
available RAMSES code Teyssier (2002), 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) DESI Collaboration
et al. (2023), 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) Grainge et al. (2017).

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 Abbott et al. (2016). 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 Corda (2009).

This research received no external funding.

Not applicable.

Not applicable.

No new data were created or analysed in this study. Data sharing is not
applicable to this article.

###### Acknowledgements.

I sincerely thank the three referees whose suggestions and corrections
improved my review. Part of this work is the result of my Ph.D. activity at
the Physics Department of the University of Turin (2017–2021). I sincerely
thank all the people that helped me during my Ph.D., in particular, Antonaldo
Diaferio, my Ph.D. supervisor, Titos Matsakos, who formulated the RG theory
together with Antonaldo Diaferio, Garry Angus, who collaborated in the study
of RG in DMS galaxies, and Andrea Sanna, who formulated CRG, together with
Antonaldo Diaferio, and Titos Matsakos. This work has been supported by the
Spoke 1 “FutureHPC & Big-Data” of the ICSC-Centro Nazionale di Ricerca in High
Performance Computing, Big Data and Quantum Computing-and hosting entity,
funded by European Union-NextGenerationEU. The author declares no conflict of
interest. Abbreviations The following abbreviations are used in this
manuscript:
BTFR | baryonic Tully–Fisher relation
---|---
CFHTLS | Canada-France-Hawaii Telescope Legacy Survey
CMB | cosmic microwave background
CRG | covariant refracted gravity
DE | dark energy
DES | Dark Energy Survey
DESI | Dark Energy Spectroscopic Instrument
DM | dark matter
DMS | DiskMass Survey
dSph | dwarf spheroidal
FLRW | Friedmann–Lema$\hat{i}$tre–Robertson–Walker
GCs | globular clusters
GR | general relativity
HSB | high surface brightness
ICM | intracluster medium
KIDS | Kilo-Degree Survey
$\Lambda$CDM |  $\Lambda$ cold dark matter
---|---
LSB | low surface brightness
MDAR | mass discrepancy–acceleration relation
MGE | multi-Gaussian expansion
MOND | modified Newtonian dynamics
RAR | radial acceleration relation
RG | refracted gravity
RMS | root-mean-square
SKA | Square Kilometer Array
SMBH | supermassive black hole
SNe Ia | Ia supernovae
SPARC | Spitzer Photometry and Accurate Rotation Curves
SPS | stellar population synthesis
TeVeS | tensor vector scalar gravity
WFL | weak field limit
References

## References

  * Planck Collaboration et al. (2020) Planck Collaboration; Aghanim, N.; Akrami, Y.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; Bartolo, N.; et al.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys. 2020, 641, A6. [CrossRef]
  * Davis et al. (1985) Davis, M.; Efstathiou, G.; Frenk, C.S.; White, S.D.M.  The evolution of large-scale structure in a universe dominated by cold dark matter.  Astrophys. J. 1985, 292, 371–394. [CrossRef]
  * Springel et al. (2006) Springel, V.; Frenk, C.S.; White, S.D.M.  The large-scale structure of the Universe.  Nature 2006, 440, 1137–1144. [CrossRef] [PubMed]
  * Markevitch (2006) Markevitch, M.  Chandra Observation of the Most Interesting Cluster in the Universe.  In _Proceedings of the The X-ray Universe 2005_ ; Wilson, A., Ed.; ESA Special Publication: Paris, France, 2006; Volume 604, p. 723, arXiv:astro-ph/0511345. [CrossRef]
  * Clowe et al. (2006) Clowe, D.; Bradač, M.; Gonzalez, A.H.; Markevitch, M.; Randall, S.W.; Jones, C.; Zaritsky, D.  A Direct Empirical Proof of the Existence of Dark Matter.  Astrophys. J. 2006, 648, L109–L113. [CrossRef]
  * Paraficz et al. (2016) Paraficz, D.; Kneib, J.P.; Richard, J.; Morandi, A.; Limousin, M.; Jullo, E.; Martinez, J.  The Bullet cluster at its best: Weighing stars, gas, and dark matter.  Astron. Astrophys. 2016, 594, A121. [CrossRef]
  * Cyburt et al. (2016) Cyburt, R.H.; Fields, B.D.; Olive, K.A.; Yeh, T.H.  Big bang nucleosynthesis: Present status.  Rev. Mod. Phys. 2016, 88, 015004. [CrossRef]
  * Zwicky (1933) Zwicky, F.  Die Rotverschiebung von extragalaktischen Nebeln.  Helv. Phys. Acta 1933, 6, 110–127.
  * Rubin and Ford (1970) Rubin, V.C.; Ford, W.K., Jr.  Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions.  Astrophys. J. 1970, 159, 379. [CrossRef]
  * Bosma (1978) Bosma, A.  The Distribution and Kinematics of Neutral Hydrogen in Spiral Galaxies of Various Morphological Types.  Ph.D. Thesis, University of Groningen, Groningen, The Netherlands, 1978.
  * Sanders (1990) Sanders, R.H.  Mass discrepancies in galaxies: Dark matter and alternatives.  Astron. Astrophys. Rev. 1990, 2, 1–28. [CrossRef]
  * Dodelson et al. (1996) Dodelson, S.; Gates, E.I.; Turner, M.S.  Cold Dark Matter.  Science 1996, 274, 69–75. [CrossRef]
  * Kirshner (1996) Kirshner, R.  Measuring the Universe with Supernovae.  In Proceedings of the APS Meeting Abstracts, Indianapolis, IN, USA, 2–5 May 1996; APS/AAPT Joint Meeting; p. F4.01.
  * Peebles and Ratra (2003) Peebles, P.J.; Ratra, B.  The cosmological constant and dark energy.  Rev. Mod. Phys. 2003, 75, 559–606. [CrossRef]
  * Weinberg (1989) Weinberg, S.  The cosmological constant problem.  Rev. Mod. Phys. 1989, 61, 1–23. [CrossRef]
  * Luo (2014) Luo, M.J.  The cosmological constant problem and re-interpretation of time.  Nucl. Phys. B 2014, 884, 344–356. [CrossRef]
  * Velten et al. (2014) Velten, H.E.S.; vom Marttens, R.F.; Zimdahl, W.  Aspects of the cosmological “coincidence problem”.  Eur. Phys. J. C 2014, 74, 3160. [CrossRef]
  * Fleury et al. (2013) Fleury, P.; Dupuy, H.; Uzan, J.P.  Can All Cosmological Observations Be Accurately Interpreted with a Unique Geometry?  Phys. Rev. Lett. 2013, 111, 091302. [CrossRef] [PubMed]
  * Douspis et al. (2019) Douspis, M.; Salvati, L.; Aghanim, N.  On the tension between Large Scale Structures and Cosmic Microwave Background.  arXiv 2019, arXiv:1901.05289.
  * Del Popolo and Le Delliou (2017) Del Popolo, A.; Le Delliou, M.  Small Scale Problems of the $\Lambda$CDM Model: A Short Review.  Galaxies 2017, 5, 17. [CrossRef]
  * de Martino et al. (2020) de Martino, I.; Chakrabarty, S.S.; Cesare, V.; Gallo, A.; Ostorero, L.; Diaferio, A.  Dark Matters on the Scale of Galaxies.  Universe 2020, 6, 107. [CrossRef]
  * van Albada and Sancisi (1986) van Albada, T.S.; Sancisi, R.  Dark Matter in Spiral Galaxies.  Philos. Trans. R. Soc. Lond. Ser. A 1986, 320, 447–464. [CrossRef]
  * Sackett (1997) Sackett, P.D.  Does the Milky Way Have a Maximal Disk?  Astrophys. J. 1997, 483, 103–110. [CrossRef]
  * Courteau and Rix (1998) Courteau, S.; Rix, H.W.  Maximal Disks and the Tully-Fisher Relation.  In _Proceedings of the Galactic Halos_ ; Zaritsky, D., Ed.; Astronomical Society of the Pacific Conference Series; Society of the Pacific Conference: San Francisco, CA, USA, 1998; Volume 136, p. 196, arXiv:astro-ph/9707290. [CrossRef]
  * Bissantz and Gerhard (2002) Bissantz, N.; Gerhard, O.  Spiral arms, bar shape and bulge microlensing in the Milky Way.  Mon. Not. R. Astron. Soc. 2002, 330, 591–608. [CrossRef]
  * Sellwood and Debattista (2014) Sellwood, J.A.; Debattista, V.P.  Re-interpretation of “Bar slowdown and the distribution of dark matter in barred galaxies” by Athanassoula.  arXiv 2014, arXiv:1410.0834. [CrossRef]
  * McGaugh and Schombert (2015) McGaugh, S.S.; Schombert, J.M.  Weighing Galaxy Disks with the Baryonic Tully-Fisher Relation.  Astrophys. J. 2015, 802, 18. [CrossRef]
  * Dubinski and Carlberg (1991) Dubinski, J.; Carlberg, R.G.  The Structure of Cold Dark Matter Halos.  Astrophys. J. 1991, 378, 496. [CrossRef]
  * Navarro et al. (1996) Navarro, J.F.; Frenk, C.S.; White, S.D.M.  The Structure of Cold Dark Matter Halos.  Astrophys. J. 1996, 462, 563. [CrossRef]
  * Mateo (1998) Mateo, M.L.  Dwarf Galaxies of the Local Group.  ARA&A 1998, 36, 435–506. [CrossRef]
  * Torrealba et al. (2019) Torrealba, G.; Belokurov, V.; Koposov, S.E.; Li, T.S.; Walker, M.G.; Sanders, J.L.; Geringer-Sameth, A.; Zucker, D.B.; Kuehn, K.; Evans, N.W.; et al.  The hidden giant: Discovery of an enormous Galactic dwarf satellite in Gaia DR2.  Mon. Not. R. Astron. Soc. 2019, 488, 2743–2766. [CrossRef]
  * Baumgardt et al. (2005) Baumgardt, H.; Grebel, E.K.; Kroupa, P.  Using distant globular clusters as a test for gravitational theories.  Mon. Not. R. Astron. Soc. 2005, 359, L1–L4. [CrossRef]
  * Baumgardt et al. (2009) Baumgardt, H.; Côté, P.; Hilker, M.; Rejkuba, M.; Mieske, S.; Djorgovski, S.G.; Stetson, P.  The velocity dispersion and mass-to-light ratio of the remote halo globular cluster NGC2419.  Mon. Not. R. Astron. Soc. 2009, 396, 2051–2060. [CrossRef]
  * Jordi et al. (2009) Jordi, K.; Grebel, E.K.; Hilker, M.; Baumgardt, H.; Frank, M.; Kroupa, P.; Haghi, H.; Côté, P.; Djorgovski, S.G.  Testing Fundamental Physics with Distant Star Clusters: Analysis of Observational Data on Palomar 14.  Astron. J. 2009, 137, 4586–4596. [CrossRef]
  * Sollima and Nipoti (2010) Sollima, A.; Nipoti, C.  Globular clusters in modified Newtonian dynamics: Velocity dispersion profiles from self-consistent models.  Mon. Not. R. Astron. Soc. 2010, 401, 131–142. [CrossRef]
  * Ibata et al. (2011a) Ibata, R.; Sollima, A.; Nipoti, C.; Bellazzini, M.; Chapman, S.C.; Dalessandro, E.  The Globular Cluster NGC 2419: A Crucible for Theories of Gravity.  Astrophys. J. 2011, 738, 186. [CrossRef]
  * Ibata et al. (2011b) Ibata, R.; Sollima, A.; Nipoti, C.; Bellazzini, M.; Chapman, S.C.; Dalessandro, E.  Polytropic Model Fits to the Globular Cluster NGC 2419 in Modified Newtonian Dynamics.  Astrophys. J. 2011, 743, 43. [CrossRef]
  * Frank et al. (2012) Frank, M.J.; Hilker, M.; Baumgardt, H.; Côté, P.; Grebel, E.K.; Haghi, H.; Küpper, A.H.W.; Djorgovski, S.G.  The velocity dispersion and mass function of the outer halo globular cluster Palomar 4.  Mon. Not. R. Astron. Soc. 2012, 423, 2917–2932. [CrossRef]
  * McGaugh et al. (2000) McGaugh, S.S.; Schombert, J.M.; Bothun, G.D.; de Blok, W.J.G.  The Baryonic Tully-Fisher Relation.  Astophys. J. 2000, 533, L99–L102. [CrossRef]
  * McGaugh (2004) McGaugh, S.S.  The Mass Discrepancy-Acceleration Relation: Disk Mass and the Dark Matter Distribution.  Astrophys. J. 2004, 609, 652–666. [CrossRef]
  * McGaugh et al. (2016) McGaugh, S.S.; Lelli, F.; Schombert, J.M.  Radial Acceleration Relation in Rotationally Supported Galaxies.  Phys. Rev. Lett. 2016, 117, 201101. [CrossRef]
  * Milgrom (1983a) Milgrom, M.  A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis.  Astrophys. J. 1983, 270, 365–370. [CrossRef]
  * Milgrom (1983b) Milgrom, M.  A modification of the Newtonian dynamics—Implications for galaxies.  Astrophys. J. 1983, 270, 371–383. [CrossRef]
  * Milgrom (1983c) Milgrom, M.  A modification of the newtonian dynamics: Implications for galaxy systems.  Astrophys. J. 1983, 270, 384–389. [CrossRef]
  * Matsakos and Diaferio (2016) Matsakos, T.; Diaferio, A.  Dynamics of galaxies and clusters in refracted gravity.  arXiv 2016, arXiv:1603.04943.
  * Sanna et al. (2023) Sanna, A.P.; Matsakos, T.; Diaferio, A.  Covariant formulation of refracted gravity.  Astron. Astrophys. 2023, 674, A209. [CrossRef]
  * Cesare et al. (2022) Cesare, V.; Diaferio, A.; Matsakos, T.  The dynamics of three nearby E0 galaxies in refracted gravity.  Astron. Astrophys. 2022, 657, A133. [CrossRef]
  * Cesare et al. (2020) Cesare, V.; Diaferio, A.; Matsakos, T.; Angus, G.  Dynamics of DiskMass Survey galaxies in refracted gravity.  Astron. Astrophys. 2020, 637, A70. [CrossRef]
  * Deur (2014) Deur, A.  A relation between the dark mass of elliptical galaxies and their shape.  Mon. Not. R. Astron. Soc. 2014, 438, 1535–1551. [CrossRef]
  * Deur (2020) Deur, A.  A correlation between the dark content of elliptical galaxies and their ellipticity.  arXiv 2020, arXiv:2010.06692.
  * McGaugh (2012) McGaugh, S.S.  The Baryonic Tully-Fisher Relation of Gas-rich Galaxies as a Test of $\Lambda$CDM and MOND.  Astron. J. 2012, 143, 40. [CrossRef]
  * Famaey and McGaugh (2012) Famaey, B.; McGaugh, S.S.  Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions.  Living Rev. Relativ. 2012, 15, 10. [CrossRef] [PubMed]
  * Bershady et al. (2010) Bershady, M.A.; Verheijen, M.A.W.; Swaters, R.A.; Andersen, D.R.; Westfall, K.B.; Martinsson, T.  The DiskMass Survey. I. Overview.  Astrophys. J. 2010, 716, 198–233. [CrossRef]
  * Cesare (2021) Cesare, V.  Dynamics of Disk and Elliptical Galaxies in Refracted Gravity.  Phys. Sci. Forum 2021, 2, 34. [CrossRef]
  * Sancisi (2004) Sancisi, R.  The visible matter—Dark matter coupling.  In Proceedings of the Dark Matter in Galaxies, Sydney, Australia, 21–25 July 2003; Ryder, S., Pisano, D., Walker, M., Freeman, K., Eds.; Astronomical Society of the Pacific: San Francisco, CA, USA, 2004; Volume 220, p. 233, arXiv:astro-ph/0311348. [CrossRef]
  * de Vaucouleurs (1959) de Vaucouleurs, G.  General Physical Properties of External Galaxies.  Handb. Der Phys. 1959, 53, 311. [CrossRef]
  * Freeman (1970) Freeman, K.C.  On the Disks of Spiral and S0 Galaxies.  Astrophys. J. 1970, 160, 811. [CrossRef]
  * Sérsic (1963) Sérsic, J.L.  Influence of the atmospheric and instrumental dispersion on the brightness distribution in a galaxy.  Bol. Asoc. Argent. Astron. Plata Argent. 1963, 6, 41–43.
  * Young (1954) Young, D.  Iterative Methods for Solving Partial Difference Equations of Elliptic Type.  Trans. Am. Math. Soc. 1954, 76, 92–111. [CrossRef]
  * Bell and de Jong (2001) Bell, E.F.; de Jong, R.S.  Stellar Mass-to-Light Ratios and the Tully-Fisher Relation.  Astrophys. J. 2001, 550, 212–229. [CrossRef]
  * Bershady et al. (2010) Bershady, M.A.; Verheijen, M.A.W.; Westfall, K.B.; Andersen, D.R.; Swaters, R.A.; Martinsson, T.  The DiskMass Survey. II. Error Budget.  Astrophys. J. 2010, 716, 234–268. [CrossRef]
  * Angus et al. (2015) Angus, G.W.; Gentile, G.; Swaters, R.; Famaey, B.; Diaferio, A.; McGaugh, S.S.; Heyden, K.J.v.d.  Mass models of disc galaxies from the DiskMass Survey in modified Newtonian dynamics.  Mon. Not. R. Astron. Soc. 2015, 451, 3551–3580. [CrossRef]
  * Milgrom (2010) Milgrom, M.  Quasi-linear formulation of MOND.  Mon. Not. R. Astron. Soc. 2010, 403, 886–895. [CrossRef]
  * Aniyan et al. (2016) Aniyan, S.; Freeman, K.C.; Gerhard, O.E.; Arnaboldi, M.; Flynn, C.  The influence of a kinematically cold young component on disc-halo decompositions in spiral galaxies: Insights from solar neighbourhood K-giants.  Mon. Not. R. Astron. Soc. 2016, 456, 1484–1494. [CrossRef]
  * Cesare et al. (2020) Cesare, V.; Colonnelli, I.; Aldinucci, M.  Practical Parallelization of Scientific Applications.  In Proceedings of the 2020 28th Euromicro International Conference on Parallel, Distributed and Network-Based Processing (PDP), Vasteras, Sweden, 11–13 March 2020; pp. 376–384. [CrossRef]
  * Aldinucci et al. (2021) Aldinucci, M.; Cesare, V.; Colonnelli, I.; Martinelli, A.R.; Mittone, G.; Cantalupo, B.; Cavazzoni, C.; Drocco, M.  Practical parallelization of scientific applications with OpenMP, OpenACC and MPI.  J. Parallel Distrib. Comput. 2021, 157, 13–29. [CrossRef]
  * Lelli et al. (2016) Lelli, F.; McGaugh, S.S.; Schombert, J.M.  SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves.  Astron. J. 2016, 152, 157. [CrossRef]
  * Li et al. (2018) Li, P.; Lelli, F.; McGaugh, S.; Schombert, J.  Fitting the radial acceleration relation to individual SPARC galaxies.  Astron. Astrophys. 2018, 615, A3. [CrossRef]
  * Lelli et al. (2017) Lelli, F.; McGaugh, S.S.; Schombert, J.M.; Pawlowski, M.S.  One Law to Rule Them All: The Radial Acceleration Relation of Galaxies.  Astrophys. J. 2017, 836, 152. [CrossRef]
  * Pota et al. (2013) Pota, V.; Forbes, D.A.; Romanowsky, A.J.; Brodie, J.P.; Spitler, L.R.; Strader, J.; Foster, C.; Arnold, J.A.; Benson, A.; Blom, C.; et al.  The SLUGGS Survey: Kinematics for over 2500 globular clusters in 12 early-type galaxies.  Mon. Not. R. Astron. Soc. 2013, 428, 389–420. [CrossRef]
  * Brodie et al. (2014) Brodie, J.P.; Romanowsky, A.J.; Strader, J.; Forbes, D.A.; Foster, C.; Jennings, Z.G.; Pastorello, N.; Pota, V.; Usher, C.; Blom, C.; et al.  The SAGES Legacy Unifying Globulars and GalaxieS Survey (SLUGGS): Sample Definition, Methods, and Initial Results.  Astrophys. J. 2014, 796, 52. [CrossRef]
  * Forbes et al. (2017) Forbes, D.A.; Alabi, A.; Brodie, J.P.; Romanowsky, A.J.; Strader, J.; Foster, C.; Usher, C.; Spitler, L.; Bellstedt, S.; Pastorello, N.; et al.  The SLUGGS Survey: A Catalog of Over 4000 Globular Cluster Radial Velocities in 27 Nearby Early-type Galaxies.  Astron. J. 2017, 153, 114. [CrossRef]
  * Mathews and Brighenti (2003) Mathews, W.G.; Brighenti, F.  Hot Gas in and around Elliptical Galaxies.  ARA&A 2003, 41, 191–239. [CrossRef]
  * Pulsoni et al. (2018) Pulsoni, C.; Gerhard, O.; Arnaboldi, M.; Coccato, L.; Longobardi, A.; Napolitano, N.R.; Moylan, E.; Narayan, C.; Gupta, V.; Burkert, A.; et al.  The extended Planetary Nebula Spectrograph (ePN.S) early-type galaxy survey: The kinematic diversity of stellar halos and the relation between halo transition scale and stellar mass.  Astron. Astrophys. 2018, 618, A94. [CrossRef]
  * Pota et al. (2015) Pota, V.; Romanowsky, A.J.; Brodie, J.P.; Peñarrubia, J.; Forbes, D.A.; Napolitano, N.R.; Foster, C.; Walker, M.G.; Strader, J.; Roediger, J.C.  The SLUGGS survey: Multipopulation dynamical modelling of the elliptical galaxy NGC 1407 from stars and globular clusters.  Mon. Not. R. Astron. Soc. 2015, 450, 3345–3358. [CrossRef]
  * Scott et al. (2013) Scott, N.; Cappellari, M.; Davies, R.L.; Verdoes Kleijn, G.; Bois, M.; Alatalo, K.; Blitz, L.; Bournaud, F.; Bureau, M.; Crocker, A.; et al.  The ATLAS3D project—XXI. Correlations between gradients of local escape velocity and stellar populations in early-type galaxies.  Mon. Not. R. Astron. Soc. 2013, 432, 1894–1913. [CrossRef]
  * Emsellem et al. (1994) Emsellem, E.; Monnet, G.; Bacon, R.  The multi-gaussian expansion method: A tool for building realistic photometric and kinematical models of stellar systems I. The formalism.  Astron. Astrophys. 1994, 285, 723–738.
  * Strader et al. (2011) Strader, J.; Romanowsky, A.J.; Brodie, J.P.; Spitler, L.R.; Beasley, M.A.; Arnold, J.A.; Tamura, N.; Sharples, R.M.; Arimoto, N.  Wide-field Precision Kinematics of the M87 Globular Cluster System.  Astrophys. J. Suppl. Ser. 2011, 197, 33. [CrossRef]
  * Zhang et al. (2007) Zhang, Z.; Xu, H.; Wang, Y.; An, T.; Xu, Y.; Wu, X.P.  Probing the Mass Distributions in NGC 1407 and Its Associated Group with the X-Ray Imaging Spectroscopic and Optical Photometric and Line-Strength Indices Data.  Astrophys. J. 2007, 656, 805–817. [CrossRef]
  * Paggi et al. (2017) Paggi, A.; Kim, D.W.; Anderson, C.; Burke, D.; D’Abrusco, R.; Fabbiano, G.; Fruscione, A.; Gokas, T.; Lauer, J.; McCollough, M.; et al.  Constraining the Physical State of the Hot Gas Halos in NGC 4649 and NGC 5846.  Astrophys. J. 2017, 844, 5. [CrossRef]
  * Fabricant et al. (1980) Fabricant, D.; Lecar, M.; Gorenstein, P.  X-ray measurements of the mass of M 87.  Astrophys. J. 1980, 241, 552–560. [CrossRef]
  * Rusli et al. (2013) Rusli, S.P.; Erwin, P.; Saglia, R.P.; Thomas, J.; Fabricius, M.; Bender, R.; Nowak, N.  Depleted Galaxy Cores and Dynamical Black Hole Masses.  Astron. J. 2013, 146, 160. [CrossRef]
  * Jeans (1915) Jeans, J.H.  On the theory of star-streaming and the structure of the universe.  Mon. Not. R. Astron. Soc. 1915, 76, 70–84. [CrossRef]
  * Mamon and Łokas (2005) Mamon, G.A.; Łokas, E.L.  Dark matter in elliptical galaxies—II. Estimating the mass within the virial radius.  Mon. Not. R. Astron. Soc. 2005, 363, 705–722. [CrossRef]
  * Cappellari (2008) Cappellari, M.  Measuring the inclination and mass-to-light ratio of axisymmetric galaxies via anisotropic Jeans models of stellar kinematics.  Mon. Not. R. Astron. Soc. 2008, 390, 71–86. [CrossRef]
  * Thomas et al. (2014) Thomas, J.; Saglia, R.P.; Bender, R.; Erwin, P.; Fabricius, M.  The Dynamical Fingerprint of Core Scouring in Massive Elliptical Galaxies.  Astrophys. J. 2014, 782, 39. [CrossRef]
  * Rantala et al. (2019) Rantala, A.; Johansson, P.H.; Naab, T.; Thomas, J.; Frigo, M.  The Simultaneous Formation of Cored, Tangentially Biased, and Kinematically Decoupled Centers in Massive Early-type Galaxies.  Astophys. J. Lett. 2019, 872, L17. [CrossRef]
  * Humphrey et al. (2006) Humphrey, P.J.; Buote, D.A.; Gastaldello, F.; Zappacosta, L.; Bullock, J.S.; Brighenti, F.; Mathews, W.G.  A Chandra View of Dark Matter in Early-Type Galaxies.  Astrophys. J. 2006, 646, 899–918. [CrossRef]
  * Bell et al. (2003) Bell, E.F.; McIntosh, D.H.; Katz, N.; Weinberg, M.D.  The Optical and Near-Infrared Properties of Galaxies. I. Luminosity and Stellar Mass Functions.  Astrophys. J. Suppl. Ser. 2003, 149, 289–312. [CrossRef]
  * Zibetti et al. (2009) Zibetti, S.; Charlot, S.; Rix, H.W.  Resolved stellar mass maps of galaxies—I. Method and implications for global mass estimates.  Mon. Not. R. Astron. Soc. 2009, 400, 1181–1198. [CrossRef]
  * Cappellari et al. (2011) Cappellari, M.; Emsellem, E.; Krajnovic, D.; McDermid, R.M.; Scott, N.; Verdoes Kleijn, G.A.; Young, L.M.; Alatalo, K.; Bacon, R.; Blitz, L.; et al.  The ATLAS3D project—I. A volume-limited sample of 260 nearby early-type galaxies: Science goals and selection criteria.  Mon. Not. R. Astron. Soc. 2011, 413, 813–836. [CrossRef]
  * Vikhlinin et al. (2006) Vikhlinin, A.; Kravtsov, A.; Forman, W.; Jones, C.; Markevitch, M.; Murray, S.S.; Van Speybroeck, L.  Chandra Sample of Nearby Relaxed Galaxy Clusters: Mass, Gas Fraction, and Mass-Temperature Relation.  Astrophys. J. 2006, 640, 691–709. [CrossRef]
  * McGaugh (2005) McGaugh, S.S.  The Baryonic Tully-Fisher Relation of Galaxies with Extended Rotation Curves and the Stellar Mass of Rotating Galaxies.  Astrophys. J. 2005, 632, 859–871. [CrossRef]
  * Sanders (2003) Sanders, R.H.  Clusters of galaxies with modified Newtonian dynamics.  Mon. Not. R. Astron. Soc. 2003, 342, 901–908. [CrossRef]
  * McGaugh (2015) McGaugh, S.S.  A tale of two paradigms: The mutual incommensurability of $\Lambda$CDM and MOND.  Can. J. Phys. 2015, 93, 250–259. [CrossRef]
  * Bekenstein (2004) Bekenstein, J.D.  Relativistic gravitation theory for the modified Newtonian dynamics paradigm.  Phys. Rev. D 2004, 70, 083509. [CrossRef]
  * Skordis (2009) Skordis, C.  TOPICAL REVIEW: The tensor-vector-scalar theory and its cosmology.  Class. Quantum Gravity 2009, 26, 143001. [CrossRef]
  * Bekenstein (2011) Bekenstein, J.D.  Tensor-vector-scalar-modified gravity: From small scale to cosmology.  Philos. Trans. R. Soc. Lond. Ser. A 2011, 369, 5003–5017. [CrossRef] [PubMed]
  * Hernandez et al. (2019) Hernandez, X.; Sussman, R.A.; Nasser, L.  Approaching the Dark Sector through a bounding curvature criterion.  Mon. Not. R. Astron. Soc. 2019, 483, 147–151. [CrossRef]
  * Skordis and Złośnik (2021) Skordis, C.; Złośnik, T.  New Relativistic Theory for Modified Newtonian Dynamics.  Phys. Rev. Lett. 2021, 127, 161302. [CrossRef] [PubMed]
  * Kunz et al. (2009) Kunz, M.; Liddle, A.R.; Parkinson, D.; Gao, C.  Constraining the dark fluid.  Phys. Rev. D 2009, 80, 083533. [CrossRef]
  * Kunz (2009) Kunz, M.  Degeneracy between the dark components resulting from the fact that gravity only measures the total energy-momentum tensor.  Phys. Rev. D 2009, 80, 123001. [CrossRef]
  * Makler et al. (2003) Makler, M.; de Oliveira, S.Q.; Waga, I.  Constraints on the generalized Chaplygin gas from supernovae observations.  Phys. Lett. B 2003, 555, 1–6. [CrossRef]
  * Bento et al. (2003) Bento, M.C.; Bertolami, O.; Sen, A.A.  Letter: Generalized Chaplygin Gas Model: Dark Energy-Dark Matter Unification and CMBR Constraints.  Gen. Relativ. Gravit. 2003, 35, 2063–2069. [CrossRef]
  * Scherrer (2004) Scherrer, R.J.  Purely Kinetic k Essence as Unified Dark Matter.  Phys. Rev. Lett. 2004, 93, 011301. [CrossRef]
  * Bruni et al. (2013) Bruni, M.; Lazkoz, R.; Rozas-Fernández, A.  Phenomenological models for unified dark matter with fast transition.  Mon. Not. R. Astron. Soc. 2013, 431, 2907–2916. [CrossRef]
  * Cadoni et al. (2018) Cadoni, M.; Casadio, R.; Giusti, A.; Mück, W.; Tuveri, M.  Effective fluid description of the dark universe.  Phys. Lett. B 2018, 776, 242–248. [CrossRef]
  * Sotiriou and Faraoni (2010) Sotiriou, T.P.; Faraoni, V.  f(R) theories of gravity.  Rev. Mod. Phys. 2010, 82, 451–497. [CrossRef]
  * Sebastiani et al. (2016) Sebastiani, L.; Vagnozzi, S.; Myrzakulov, R.  Mimetic gravity: A review of recent developments and applications to cosmology and astrophysics.  arXiv 2016, arXiv:1612.08661. [CrossRef]
  * Berezhiani et al. (2017) Berezhiani, L.; Khoury, J.; Wang, J.  Universe without dark energy: Cosmic acceleration from dark matter-baryon interactions.  Phys. Rev. D 2017, 95, 123530. [CrossRef]
  * Ferreira et al. (2019) Ferreira, E.G.M.; Franzmann, G.; Khoury, J.; Brandenberger, R.  Unified superfluid dark sector.  J. Cosmol. Astropart. Phys. 2019, 2019, 027. [CrossRef]
  * Arbey and Coupechoux (2021) Arbey, A.; Coupechoux, J.F.  Unifying dark matter, dark energy and inflation with a fuzzy dark fluid.  J. Cosmol. Astropart. Phys. 2021, 2021, 033. [CrossRef]
  * Campigotto et al. (2019) Campigotto, M.C.; Diaferio, A.; Fatibene, L.  Conformal gravity: Light deflection revisited and the galactic rotation curve failure.  Class. Quantum Gravity 2019, 36, 245014. [CrossRef]
  * Suzuki et al. (2012) Suzuki, N.; Rubin, D.; Lidman, C.; Aldering, G.; Amanullah, R.; Barbary, K.; Barrientos, L.F.; Botyanszki, J.; Brodwin, M.; Connolly, N.; et al.  The Hubble Space Telescope Cluster Supernova Survey. V. Improving the Dark-energy Constraints above z > 1 and Building an Early-type-hosted Supernova Sample.  Astrophys. J. 2012, 746, 85. [CrossRef]
  * Amendola and Tsujikawa (2015) Amendola, L.; Tsujikawa, S.  Dark Energy: Theory and Observations; Cambridge University Press: Cambridge, UK, 2015.
  * Copeland et al. (2006) Copeland, E.J.; Sami, M.; Tsujikawa, S.  Dynamics of Dark Energy.  Int. J. Mod. Phys. D 2006, 15, 1753–1935. [CrossRef]
  * Frusciante and Perenon (2020) Frusciante, N.; Perenon, L.  Effective field theory of dark energy: A review.  Phys. Rep. 2020, 857, 1–63. [CrossRef]
  * Wen et al. (2018) Wen, S.; Wang, S.; Luo, X.  Comparing dark energy models with current observational data.  J. Cosmol. Astropart. Phys. 2018, 2018, 011. [CrossRef]
  * Capozziello et al. (2006) Capozziello, S.; Nojiri, S.; Odintsov, S.D.  Dark energy: The equation of state description versus scalar-tensor or modified gravity.  Phys. Lett. B 2006, 634, 93–100. [CrossRef]
  * Gerardi et al. (2019) Gerardi, F.; Martinelli, M.; Silvestri, A.  Reconstruction of the Dark Energy equation of state from latest data: The impact of theoretical priors.  J. Cosmol. Astropart. Phys. 2019, 2019, 042. [CrossRef]
  * Chevallier and Polarski (2001) Chevallier, M.; Polarski, D.  Accelerating Universes with Scaling Dark Matter.  Int. J. Mod. Phys. D 2001, 10, 213–223. [CrossRef]
  * Linder (2003) Linder, E.V.  Exploring the Expansion History of the Universe.  Phys. Rev. Lett. 2003, 90, 091301. [CrossRef] [PubMed]
  * Bean and Melchiorri (2002) Bean, R.; Melchiorri, A.  Current constraints on the dark energy equation of state.  Phys. Rev. D 2002, 65, 041302. [CrossRef]
  * Solà and Štefančić (2005) Solà, J.; Štefančić, H.  Effective equation of state for dark energy: Mimicking quintessence and phantom energy through a variable $\Lambda$.  Phys. Lett. B 2005, 624, 147–157. [CrossRef]
  * Hazra et al. (2015) Hazra, D.K.; Majumdar, S.; Pal, S.; Panda, S.; Sen, A.A.  Post-Planck dark energy constraints.  Phys. Rev. D 2015, 91, 083005. [CrossRef]
  * Di Paolo et al. (2019) Di Paolo, C.; Salucci, P.; Fontaine, J.P.  The Radial Acceleration Relation (RAR): Crucial Cases of Dwarf Disks and Low-surface-brightness Galaxies.  Astrophys. J. 2019, 873, 106. [CrossRef]
  * Santos-Santos et al. (2020) Santos-Santos, I.M.E.; Navarro, J.F.; Robertson, A.; Benítez-Llambay, A.; Oman, K.A.; Lovell, M.R.; Frenk, C.S.; Ludlow, A.D.; Fattahi, A.; Ritz, A.  Baryonic clues to the puzzling diversity of dwarf galaxy rotation curves.  Mon. Not. R. Astron. Soc. 2020, 495, 58–77. [CrossRef]
  * Jiménez et al. (2013) Jiménez, M.A.; Garcia, G.; Hernandez, X.; Nasser, L.  The Massive Elliptical Galaxy NGC 4649 from the Perspective of Extended Gravity.  Astrophys. J. 2013, 768, 142. [CrossRef]
  * Durazo et al. (2017) Durazo, R.; Hernandez, X.; Cervantes Sodi, B.; Sánchez, S.F.  A Universal Velocity Dispersion Profile for Pressure Supported Systems: Evidence for MONDian Gravity across Seven Orders of Magnitude in Mass.  Astrophys. J. 2017, 837, 179. [CrossRef]
  * Durazo et al. (2018) Durazo, R.; Hernandez, X.; Cervantes Sodi, B.; Sanchez, S.F.  A Test of MONDian Gravity in $\sim$300 Pressure-supported Elliptical Galaxies from the MaNGA Survey.  Astrophys. J. 2018, 863, 107. [CrossRef]
  * Salucci et al. (2012) Salucci, P.; Wilkinson, M.I.; Walker, M.G.; Gilmore, G.F.; Grebel, E.K.; Koch, A.; Frigerio Martins, C.; Wyse, R.F.G.  Dwarf spheroidal galaxy kinematics and spiral galaxy scaling laws.  Mon. Not. R. Astron. Soc. 2012, 420, 2034–2041. [CrossRef]
  * Oh et al. (2015) Oh, S.H.; Hunter, D.A.; Brinks, E.; Elmegreen, B.G.; Schruba, A.; Walter, F.; Rupen, M.P.; Young, L.M.; Simpson, C.E.; Johnson, M.C.; et al.  High-resolution Mass Models of Dwarf Galaxies from LITTLE THINGS.  Astron. J. 2015, 149, 180. [CrossRef]
  * Walker et al. (2007) Walker, M.G.; Mateo, M.; Olszewski, E.W.; Gnedin, O.Y.; Wang, X.; Sen, B.; Woodroofe, M.  Velocity Dispersion Profiles of Seven Dwarf Spheroidal Galaxies.  Astophys. J. 2007, 667, L53–L56. [CrossRef]
  * Salomon et al. (2015) Salomon, J.B.; Ibata, R.A.; Martin, N.F.; Famaey, B.  The intrinsic ellipticity of dwarf spheroidal galaxies: Constraints from the Andromeda system.  Mon. Not. R. Astron. Soc. 2015, 450, 1409–1419. [CrossRef]
  * Rines and Diaferio (2006) Rines, K.; Diaferio, A.  CIRS: Cluster Infall Regions in the Sloan Digital Sky Survey. I. Infall Patterns and Mass Profiles.  Astron. J. 2006, 132, 1275–1297. [CrossRef]
  * Rines et al. (2013) Rines, K.; Geller, M.J.; Diaferio, A.; Kurtz, M.J.  Measuring the Ultimate Halo Mass of Galaxy Clusters: Redshifts and Mass Profiles from the Hectospec Cluster Survey (HeCS).  Astrophys. J. 2013, 767, 15. [CrossRef]
  * Gianfagna et al. (2021) Gianfagna, G.; De Petris, M.; Yepes, G.; De Luca, F.; Sembolini, F.; Cui, W.; Biffi, V.; Kéruzoré, F.; Macías-Pérez, J.; Mayet, F.; et al.  Exploring the hydrostatic mass bias in MUSIC clusters: Application to the NIKA2 mock sample.  Mon. Not. R. Astron. Soc. 2021, 502, 5115–5133. [CrossRef]
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration.  The Dark Energy Survey.  arXiv 2005, arXiv:astro-ph/0510346. [CrossRef]
  * Euclid Collaboration et al. (2022) Euclid Collaboration.; Scaramella, R.; Amiaux, J.; Mellier, Y.; Burigana, C.; Carvalho, C.S.; Cuillandre, J.C.; Da Silva, A.; Derosa, A.; Dinis, J.; et al.  Euclid preparation. I. The Euclid Wide Survey.  Astron. Astrophys. 2022, 662, A112. [CrossRef]
  * Amendola et al. (2018) Amendola, L.; Appleby, S.; Avgoustidis, A.; Bacon, D.; Baker, T.; Baldi, M.; Bartolo, N.; Blanchard, A.; Bonvin, C.; Borgani, S.; et al.  Cosmology and fundamental physics with the Euclid satellite.  Living Rev. Relativ. 2018, 21, 2. [CrossRef] [PubMed]
  * Huterer et al. (2015) Huterer, D.; Kirkby, D.; Bean, R.; Connolly, A.; Dawson, K.; Dodelson, S.; Evrard, A.; Jain, B.; Jarvis, M.; Linder, E.; et al.  Growth of cosmic structure: Probing dark energy beyond expansion.  Astropart. Phys. 2015, 63, 23–41. [CrossRef]
  * Noller and Nicola (2019) Noller, J.; Nicola, A.  Cosmological parameter constraints for Horndeski scalar-tensor gravity.  Phys. Rev. D 2019, 99, 103502. [CrossRef]
  * Planck Collaboration et al. (2020) Planck Collaboration; Aghanim, N.; Akrami, Y.; Arroja, F.; Ashdown, M.; Aumont, J.; Baccigalupi, C.; Ballardini, M.; Banday, A.J.; Barreiro, R.B.; et al.  Planck 2018 results. I. Overview and the cosmological legacy of Planck.  Astron. Astrophys. 2020, 641, A1. [CrossRef]
  * de Jong et al. (2013) de Jong, J.T.A.; Kuijken, K.; Applegate, D.; Begeman, K.; Belikov, A.; Blake, C.; Bout, J.; Boxhoorn, D.; Buddelmeijer, H.; Buddendiek, A.; et al.  The Kilo-Degree Survey.  Messenger 2013, 154, 44–46. [CrossRef]
  * Gwyn (2012) Gwyn, S.D.J.  The Canada-France-Hawaii Telescope Legacy Survey: Stacked Images and Catalogs.  Astron. J. 2012, 143, 38. [CrossRef]
  * di Porto and Amendola (2008) di Porto, C.; Amendola, L.  Observational constraints on the linear fluctuation growth rate.  Phys. Rev. D 2008, 77, 083508. [CrossRef]
  * Bueno Sánchez and Perivolaropoulos (2011) Bueno Sánchez, J.C.; Perivolaropoulos, L.  Dark energy and matter perturbations in scalar-tensor theories of gravity.  J. Phys. Conf. Ser. 2011, _283_ , 012006. [CrossRef]
  * Pace et al. (2014) Pace, F.; Moscardini, L.; Crittenden, R.; Bartelmann, M.; Pettorino, V.  A comparison of structure formation in minimally and non-minimally coupled quintessence models.  Mon. Not. R. Astron. Soc. 2014, 437, 547–561. [CrossRef]
  * Kofinas and Lima (2017) Kofinas, G.; Lima, N.A.  Dynamics of cosmological perturbations in modified Brans-Dicke cosmology with matter-scalar field interaction.  Phys. Rev. D 2017, 96, 084016. [CrossRef]
  * Teyssier (2002) Teyssier, R.  Cosmological hydrodynamics with adaptive mesh refinement. A new high resolution code called RAMSES.  Astron. Astrophys. 2002, 385, 337–364. [CrossRef]
  * DESI Collaboration et al. (2023) DESI Collaboration; Adame, A.G.; Aguilar, J.; Ahlen, S.; Alam, S.; Aldering, G.; Alexander, D.M.; Alfarsy, R.; Allende Prieto, C.; Alvarez, M.; et al.  The Early Data Release of the Dark Energy Spectroscopic Instrument.  arXiv 2023, arXiv:2306.06308. [CrossRef]
  * Grainge et al. (2017) Grainge, K.; Alachkar, B.; Amy, S.; Barbosa, D.; Bommineni, M.; Boven, P.; Braddock, R.; Davis, J.; Diwakar, P.; Francis, V.; et al.  Square Kilometre Array: The radio telescope of the XXI century.  Astron. Rep. 2017, 61, 288–296. [CrossRef]
  * Abbott et al. (2016) Abbott, B.P.; Abbott, R.; Abbott, T.D.; Abernathy, M.R.; Acernese, F.; Ackley, K.; Adams, C.; Adams, T.; Addesso, P.; Adhikari, R.X.; et al.  Observation of Gravitational Waves from a Binary Black Hole Merger.  Phys. Rev. Lett. 2016, 116, 061102. [CrossRef] [PubMed]
  * Corda (2009) Corda, C.  Interferometric Detection of Gravitational Waves: The Definitive Test for General Relativity.  Int. J. Mod. Phys. D 2009, 18, 2275–2282. [CrossRef]