Dark Coincidences: Small-Scale Solutions with Refracted Gravity and MOND

General relativity and its Newtonian weak field limit are not sufficient to explain the observed phenomenology in the Universe, from the formation of large-scale structures to the dynamics of galaxies, with the only presence of baryonic matter. The most investigated cosmological model, the $\Lambda$CDM, accounts for the majority of observations by introducing two dark components, dark energy and dark matter, which represent $\sim$95% of the mass-energy budget of the Universe. Nevertheless, the $\Lambda$CDM model faces important challenges on the scale of galaxies. For example, some very tight relations between the properties of dark and baryonic matters in disk galaxies, such as the baryonic Tully-Fisher relation (BTFR), the mass discrepancy-acceleration relation (MDAR), and the radial acceleration relation (RAR), which see the emergence of the acceleration scale $a_0 \simeq 1.2 \times 10^{-10}$ m s$^{-2}$, cannot be intuitively explained by the CDM paradigm, where cosmic structures form through a stochastic merging process. An even more outstanding coincidence is due to the fact that the acceleration scale $a_0$, emerging from galaxy dynamics, also seems to be related to the cosmological constant $\Lambda$. Another challenge is provided by dwarf galaxies, which are darker than what is expected in their innermost regions. These pieces of evidence can be more naturally explained, or sometimes even predicted, by modified theories of gravity, that do not introduce any dark fluid. I illustrate possible solutions to these problems with the modified theory of gravity MOND, which departs from Newtonian gravity for accelerations smaller than $a_0$, and with Refracted Gravity, a novel classical theory of gravity introduced in 2016, where the modification of the law of gravity is instead regulated by a density scale.


Introduction
If we assume the validity of standard gravity, that is, General Relativity (GR), and the presence of ordinary (baryonic) matter alone, we cannot explain the observations from the largest to the smallest scales of the Universe. The cosmic microwave background (CMB) radiation [1], the dynamics and the large-scale distribution of the cosmic structure [2,3], the gravitational lensing, with particular reference to the Bullet Cluster [4][5][6], the dynamics of galaxy clusters [7], and the flat trend of the rotation curves of disk galaxies [8][9][10] all show a mass discrepancy of ∼80-90%. Moreover, the Hubble diagram of Ia Supernovae [11,12] proved the expansion of the Universe to be accelerated, which is not what we expect from the attractive nature of the gravity force. This probably represents the most important open question in modern cosmology and it launches several lines of research.
With the only presence of baryonic matter, GR cannot justify ∼95% of the components of the Universe. The most investigated solution is provided by the Λ cold dark matter (ΛCDM) cosmological model, which assumes the validity of GR and introduces two dark components besides baryonic matter to explain this missing ∼95%. on the galaxy scale, such as the three mentioned scaling relations and the difference between the dynamics of HSB and LSB galaxies.
Another theory of modified gravity, more recently introduced, is refracted gravity (RG) [44], a classical theory where the modification of the law of gravity is regulated by the value of the local mass density, rather than of the acceleration. RG has shown some encouraging results in describing the dynamics of galaxies of different shapes, such as disk [45] and elliptical E0 galaxies with nearly spherical morphology [46]. In RG, the shape of the gravitational field lines depends on the morphology of the system: they are refracted towards the mid-plane of flattened systems and they remain radial for spherical systems, which might intuitively explain the different dynamics of dwarf galaxies (more flattened) and globular clusters (nearly spherical). Moreover, the relativistic formulation of RG [47] belongs to the class of scalar-tensor theories and introduces a single scalar field that explains the phenomenology of both DM and DE. The covariant RG provides a natural explanation for the a 0 ∼ Λ 1/2 relation, suggesting a unification of the two dark sectors.
In this review, we present more intuitive explanations, with the theories MOND and RG, for the three scaling relations, for the emergence of the same acceleration scale from these relations and from the DE sector, and for the different dynamics of LSB, dwarf, and dSph galaxies and of GCs. The outline of the paper develops as follows. In Section 2, we detail the BTFR, the MDAR, and the RAR (Section 2.1) and we explain the different interpretations of these relations in a ΛCDM framework (Section 2.2). Section 3 describes the formulation of MOND and RG theories (Section 3.1) and the interpretation of the three scaling relations with these theories (Section 3.2). Section 4 illustrates how the a 0 acceleration scale and the a 0 -Λ relation emerge in Newtonian, MOND, and RG theories. Section 5 describes the different dynamical properties of LSB, dwarf, and dSph galaxies and of GCs and how they can be modelled or interpreted in Newtonian, MOND, and RG gravities. Section 6 concludes the paper.

The Baryonic Scaling Relations
The mass discrepancy on the galaxy scale can be neatly quantified by three relations that tightly correlate the properties of dark and baryonic matters in galaxies: the BTFR, the MDAR, and the RAR.

Description of the Three Relations
The BTFR [20] (Figure 1) correlates the total baryonic mass, M bar , and the asymptotic value of the flat part of the rotation curve of galaxies, V f , according to the relation, e.g., [48]: where the normalisation A and the slope b are free parameters. In the BTFR, each point corresponds to one galaxy. McGaugh [48] fits the (V f ,M bar ) data points from 47 gas-rich galaxies with Equation (1), adopting different techniques. The normalisations resulting from these fits are in agreement with each other and the slopes are consistent with 4. In the BTFR, the acceleration scale a 0 emerges from its normalisation. Setting the slope to 4, McGaugh [48] found the normalisation A = (47 ± 6) M km −4 s 4 , comparable to the expression (Ga 0 ) −1 .
For mass-to-light ratios in the 3.6 µm band M/L [3.6] 0.5 M /L , the BTFR intrinsic scatter is minimised to ∼0.10 dex [49]. The residuals of the measured BTFR from the model (Equation (1)) do not correlate with galaxy properties, such as the radius or the surface brightness [50].  [51]. M b is the total galaxy baryonic mass, accounting for the stars and the gas, and V f is the mean circular velocity along the flat region of the rotation curve. The black solid line represents the linear fit to the data with the relation log(M b /M ) = s log(V f /km s −1 ) + I. σ ⊥ denotes the intrinsic scatter of the BTFR, expressed in dex. The galaxies are colour-coded according to the gas fraction F g to the total baryonic mass. The figure is re-adapted from Figure 2 in [52].
Two local, rather than global, scaling relations related to the BTFR are the MDAR and the RAR [53]. The MDAR [21] anti-correlates, at each distance R from the galaxy centre, the Newtonian acceleration generated by the baryons distribution, g bar , and the squared ratio (V/V bar ) 2 , where V and V bar are the total and the baryons-only velocities ( Figure 2, top panel). If we assume spherical symmetry, the ratio (V/V bar ) 2 coincides with the mass discrepancy, M/M bar , where M and M bar are the masses of the entire galaxy and of its baryonic component.
The MDAR can be modelled by the following relation, e.g., [54]: As long as the acceleration g bar is a 0 , the mass discrepancy maintains around 1. Instead, when g bar goes below a 0 , the mass discrepancy starts to increase. The intrinsic scatter in both the MDAR and the BTFR is minimised by the same mass-to-light ratio, consistent with the estimates from stellar population synthesis (SPS) models [21]. In fact, the MDAR can also be expressed as a relation between the mass discrepancy and g obs (Figure 2, bottom panel), that is, the total observed acceleration, which presents a slightly larger scatter than the relation with g bar [23]. This relation can be reproduced by replacing g bar with g obs in Equation (2). The mass discrepancy-acceleration relation (MDAR) represented as the squared ratio between the total (V) and the baryonic (V b ) velocities as a function of the Newtonian acceleration due to baryonic matter, g N (top panel) and the total acceleration a (bottom panel). We can clearly see that the mass discrepancy is observed from accelerations smaller than a 0 10 −10 m s −2 (indicated in the figure as a vertical blue dashed line, for reference). The (V/V b ) 2 − a relation shows a slightly larger scatter than the The black dots represent hundreds of individual resolved data points belonging to the rotation curves of about one hundred spiral galaxies. The black solid line (V/V b ) 2 = 1, showing the no mass discrepancy case, is illustrated as a reference in both panels. The figure is re-adapted from Figure 10 in [23].
A slightly different perspective with respect to the MDAR is provided by the RAR [55] ( Figure 3), which links the centripetal acceleration inferred from the observed rotation curve V(R), g obs = V 2 /R, with the Newtonian acceleration due to the baryons alone, g bar [22]. The data in the (g bar ,g obs ) plane of the RAR present an even tighter relation with respect to the data in the (g bar ,(V/V bar ) 2 ) plane of the MDAR. Moreover, the (g bar ,g obs ) plane presents an advantage compared to the (g bar ,(V/V bar ) 2 ) plane, since the g bar and g obs quantities and their corresponding uncertainties are completely independent from each other [56].
McGaugh and coauthors [22] fitted the (g bar ,g obs ) data from 153 edge-on disk galaxies belonging to Spitzer photometry and accurate rotation curves (SPARC) catalogue [51] with the relation: where the only free parameter g † = (1.20 ± 0.02 ± 0.24) × 10 −10 m s −2 is 1σ consistent with a 0 . The errors of 0.02 × 10 −10 m s −2 and of 0.24 × 10 −10 m s −2 represent the random and the systematic contributions to the uncertainty, respectively. In particular, the random error represents the 1σ confidence interval and the systematic error represents the 20% normalisation uncertainty due to the fact that the mass-to-light ratios of the disks and the bulges of the galaxies are kept fixed across the SPARC sample. Assuming a disk and a bulge mass-to-light ratio M/L [3.6] of 0.5 and 0.7 M /L , respectively, which are reasonable values in the 3.6 µm band, the RAR of SPARC galaxies is retrieved with an observed scatter of 0.13 dex [22]. This value closely coincides with the scatter of 0.12 dex due to the observational errors of the measured rotation curves, distances, and galaxy inclinations, and to the possible variation of the mass-to-light ratio among galaxies, which leaves little room for intrinsic scatter [22,53].
Li and collaborators [56] wanted to test whether the RAR was followed by the individual galaxies in the SPARC sample. They fitted Equation (3) to the observed RAR of the individual SPARC galaxies, finding a RAR with a smaller scatter of [0.054-0.057] dex, both by fixing g † to a 0 and by leaving it free to vary. Since, differently from [22], they estimated the mass-to-light ratios and marginalised over the errors on the galaxy distances and inclinations from the RAR of the single SPARC galaxies, the obtained scatter might be assimilated to the intrinsic one, rather than to the observed one. The radial acceleration relation (RAR) built from Spitzer photometry and accurate rotation curves (SPARC) galaxies. For each galaxy, a mass-to-light ratio in the 3.6 µm band for the disk and the bulge equal to 0.5 and 0.7 M /L , respectively, is adopted. The black solid line represents the fit to the data with Equation (3) and the black dashed line is the g obs = g bar relation, for reference. The only best-fit parameter, g † , is highlighted in the top-left corner of the panel and it is 1σ consistent with the acceleration scale a 0 . The histogram and the bottom panel show the distribution of the residuals of the observed RAR from Equation (3), and its standard deviation σ, which quantifies the RAR observed scatter, is shown in the top-left corner of the panel. The figure is re-adapted from Figure 2 in [57].
The three relations cover a baryonic mass range of six orders of magnitude, from M bar ∼ 10 12 M , correspondent to the most massive HSB spiral galaxies, to M bar ∼ 10 6 M , correspondent to the dwarf and the LSB galaxies. However, the observed scatter of the RAR might increase from 0.13 dex to 0.24 dex for small g bar due to dwarf and LSB galaxies with slowly-rising rotation curves [58], which requires further investigation.
The works of McGaugh et al. and of Kroupa et al. [59,60] might demonstrate a fundamental origin for the acceleration scale a 0 , and, thus, for the RAR, since they show a consistency of g † among different galaxies. The work of Li et al. [56] might suggest the same result, since fitting Equation (3) to the RAR of the individual SPARC galaxies by leaving g † free to vary does neither improve the χ 2 nor reduce the obtained scatter with respect to the fits in which g † is fixed to a 0 = 1.2 × 10 −10 m s −2 . However, several works question this result. For example, Rodrigues et al. [61] estimated the g † parameter with Bayesian inference from 193 individual galaxies from SPARC [51] and THINGS [62,63] samples, rejecting the consistency of g † among the galaxies at the 10σ level. The discrepancy reduces but does not disappear by removing some approximations, becoming equal to 5σ [64]. The work of Zhou et al. [65] leads to a similar result. If this is the case, it would represent an important challenge for MOND (see Section 3.2). A debate about the existence or the absence of a universal a 0 is presented in [66]. Another study sees the emergence of an acceleration scale consistent with a 0 , besides from the BTFR, built from rotationally-supported galaxies, from the baryonic Faber-Jackson relation built from pressure-supported galaxies, such as elliptical galaxies and GCs, and from the baryonic Faber-Jackson relation built from galaxy clusters ( [67] and references). The fact that the same acceleration scale also emerges from the dynamics of galaxy clusters, which are different from galaxies both in terms of scale and evolutionary histories, might pose an even more severe issue for ΛCDM.

Interpretation of the Three Relations in Newtonian Gravity
The ΛCDM paradigm does not provide a natural explanation for these three relations. Indeed, these relations correlate quantities related to the dynamics of a galaxy to its baryonic content, which is counterintuitive in a universe where the dynamics of structures is dominated by DM. The pieces of evidence that are most difficult to interpret are the small scatter, in agreement with observational uncertainties, and the lack of correlations of the residuals between the observed and the modelled relations with the galaxy properties. This phenomenology indicates a quite precise fine-tuning between dark and baryonic matter in galaxies. Moreover, the MDAR and the RAR pose a more serious issue for ΛCDM compared to the BTFR due to their local, rather than global, nature.
The three relations resulting from semi-analytical models of galaxy formation and DMonly cosmological simulations in ΛCDM do not completely agree with the observational data. The predicted minimum scatter of the BTFR is larger than the observed one (0.17 dex vs 0.10 dex) [68,69] and the simulated BTFR slope b is equal to 3, 8σ discrepant from the value b = 3.98 ± 0.12, fitted by McGaugh and coauthors [20] from the observed BTFR. However, some successful explanations for the three relations remaining within the ΛCDM paradigm exist, even if with some issues.
Ludlow et al. [70] demonstrated that a set of galaxies resulting from the EAGLE suite of ΛCDM hydrodynamical simulations, run with the same initial conditions but with different stellar and active galactic nuclei (AGN) feedback processes for the baryonic component, follow a RAR-like relation, independently from the considered feedback. Different feedback models make simulated galaxies moving along the RAR and not perpendicular to it, which produces a small RAR scatter 0.08 dex. However, the g † parameter fitted with Equation (3) from the galaxies simulated in [70] results in a 70σ inconsistency with a 0 = 1.2 × 10 −10 m s −2 . Moreover, measurement errors are not included in the simulated galaxies and the obtained result about the scatter should be further investigated [56]. The same consideration can be drawn for the scatter of 0.06 dex of the RAR obtained from the 32 galaxies resulting from the MUGS2 "zoomin" hydrodynamic simulations in ΛCDM [71]. On the other hand, Stone and Courteau [72] found that the intrinsic scatter of the stellar RAR, where only stellar mass is considered to compute the baryonic acceleration, is of (0.11 ± 0.02) dex, in agreement, despite being slightly larger, with ΛCDM predictions. They obtained this result from PROBES, a catalogue made of more than 2500 spiral galaxies taken from six deep imaging and spectroscopic surveys. Yet, neglecting gas masses might affect the obtained result. In fact, the question of the RAR raises further complications since for small g bar , namely for dwarf and LSB galaxies, its scatter increases [58] (see Section 2.1), and the RAR built from some galaxy samples different from SPARC shows some correlations between its residuals and certain galaxy properties [45,73].
Some hydrodynamical simulations, which are simulations that include the presence of baryons, can reproduce the slope of the BTFR [74,75]. Yet, the small scatter of the BTFR can only be explained with a quite precise balance between star formation efficiency and stellar feedback processes [48,49].
With the semi-empirical model proposed by Di Cintio and Lelli [69], the shape and the scatter of the MDAR are reproduced but this does not simultaneously account for the small scatter of ∼0.10 dex of the BTFR, which instead results in 0.17 dex.
Mayer et al. [54] used the Magneticum hydrodynamical simulation, which provides a large and representative sample of galaxies covering a large range of masses and a variety of morphologies, from z ∼ 0 to z ∼ 3, to see whether the baryonic scaling relations (BTFR, MDAR, and RAR) were reproduced in a ΛCDM context. The resulting BTFR, built for galaxies at different redshifts in the range of 0.1 < z < 2.3, has a slope more consistent with 3, rather than with 4, as in the observed BTFR. The MDAR and the RAR built from simulated Magneticum galaxies at redshift z ∼ 0.1 reproduce the observed relations (Equations (2) and (3)) with a fitted acceleration scale consistent with a 0 = 1.2 × 10 −10 m s −2 and a scatter in agreement with the observations from the SPARC sample. However, the simulated MDAR and RAR show a positive correlation between the total baryonic mass and the mass discrepancy which is not observed in SPARC data. Other ΛCDM hydrodynamical simulations of single galaxies might indicate that the baryonic scaling relations naturally emerge in this framework as well [76,77].
It is important to note that many small-scale problems of ΛCDM, such as too large bulges in disk galaxies, the cusp/core problem, and the too-big-to-fail problem, have been, at least partially, solved by introducing baryonic feedback mechanisms, such as outflows due to AGN and Supernovae, and, therefore, it is not so unlikely to suppose that galaxy dynamics is regulated by baryonic physics [54].

MOND and RG
A totally different scenario for the three scaling relations is provided by theories of modified gravity that do not imply the presence of DM. Without DM, the fine-tuning issue between the properties of dark and baryonic matters disappears. A modified gravity theory that not only describes but even predicted these scaling relations is MOND. Another theory of modified gravity that seems to account for these relations is RG. In the following subsections, I present the formulation of these two theories of gravity and how they reproduce the mentioned scaling relations.

MOND
In 1983, Milgrom [41][42][43] formulated MOND, a theory of gravity that mimics the effect of DM with a boost of the gravitational field compared to the Newtonian one in low-acceleration environments. MOND is a general paradigm that assumes spacetime scale-invariance when the acceleration is a a 0 . Specifically, the acceleration a presents the following asymptotic values: where g N is the Newtonian acceleration. The MOND paradigm is obtained by modifying either gravity, e.g., [78] or inertia [79], where the modified-inertia version [79] was less developed. In the first nonrelativistic modified gravity version of MOND [78], the following Poisson equation was defined: where the interpolating function µ, monotonic in its argument, has these two asymptotic behaviours: For a a 0 , the Newtonian Poisson equation is retrieved: and for a a 0 , we observe a boost of the gravitational field over the Newtonian one. In this regime, whereas the Newtonian field is ∝ R −2 , the MOND field is ∝ R −1 , deviating from the Newtonian inverse square law and reproducing the flat trend of the rotation curves without the presence of DM.
The MOND paradigm was also defined with other modified gravity formulations, such as QUMOND [80], where the MONDian behaviour of the gravitational field is obtained with the Poisson equation: where φ N is Newtonian gravitational potential. A possible form of the interpolating function ν is given by the "simple ν-function" (Equation (50) with n = 1 in [23]): with y being equal to

Refracted Gravity
RG is a novel theory of modified gravity inspired by the behaviour of electrodynamics in matter that does not resort to DM [44]. RG was formulated in a nonrelativistic way by [44] and its gravitational potential φ obeys the modified Poisson equation: where the gravitational permittivity (ρ) mimics the DM phenomenology. Whereas in MOND the modification of the law of gravity is regulated by an acceleration scale, in RG it is regulated by a mass density scale. The gravitational permittivity (ρ) is a monotonic increasing function of the local mass density ρ, it depends on three universal free parameters, and it presents the following asymptotic limits in the high and low-density regimes: where the permittivity of vacuum 0 and the critical density ρ c are two of the three free parameters of the theory. As in MOND for a a 0 , for ρ ρ c , Equation (10) reduces to the Newtonian Poisson equation (7). When ρ ρ c , the RG field is boosted compared to the Newtonian one thanks to the value of 0 ∈ (0, 1).
RG predicts a different behaviour for the gravitational field in the low-density environments of spherical and flattened systems. In the external regions of spherical systems, where ρ ρ c , the gravitational field does not deviate from the inverse square law, i.e., ∂φ/∂r ∝ r −2 , as in the Newtonian case, but it is boosted with respect to the Newtonian field by the inverse of the gravitational permittivity, as obtained by integrating Equation (10): where M(< r) is the system mass within the spherical radius r.
Whereas for spherical systems we still observe a Newtonian trend, in the outskirts of flattened systems the field lines are refracted toward the mid-plane of the object. This can be seen by expanding the left-hand side of Equation (10): where the term " ∂ ∂ρ ∇ρ · ∇φ", different from zero in nonspherical systems, is responsible for the focussing of the field lines. In nonspherical configurations, we thus observe an analogy with electrodynamics in matter: the gravitational field lines behave like electric field lines when they cross a dielectric medium with a nonuniform permittivity, changing both in direction and in magnitude.
This redirection effect in the low-density regions of flattened systems yields the trend ∂φ/∂R ∼ a 0 |g N | ∝ R −1 for the RG gravitational field [44], where |g N | = |∂φ N /∂R| is the Newtonian field. In this regime, the RG field deviates from the Newtonian inverse square law and is subject to a boost that in Newtonian gravity is obtained with the presence of DM. This limit also coincides with the MOND asymptotic behaviour for a a 0 (see Equation (4)) and it suggests that the ability of MOND in describing galaxy dynamics is shared by RG, as demonstrated in [45]. According to Equation (13), the flatter the system, the larger the boost of the gravitational field, and, thus, the larger the mass discrepancy, as interpreted in Newtonian theory. Figure 4 summarises the analogies and the differences between Newtonian (top panels) and RG (bottom panels) gravitational fields for flat (left panels) and spherical (right panels) systems. The behaviour of the RG gravitational field in the low-density regime is comparable to the behaviour of the MOND gravitational field in the low-acceleration regime only for nonspherical configurations. In the low-acceleration regions of spherical systems, MOND field is ∝ R −1 , as in the low-acceleration regions of flattened systems, whereas, in the low-density regions of spherical systems, the R-dependence of RG field remains Newtonian. Moreover, RG and MOND fields also present a difference for flattened systems. As we can see in the left panels of Figure 4, the refraction effect of the RG field lines in a flattened object already begins where the local density ρ is still above ρ c and, in that region, the RG field already begins to be boosted compared to Newtonian one. Instead, MOND modification with respect to the Newtonian field only appears in regions where a < a 0 .
In all the analyses of galaxy dynamics performed with RG [44][45][46], the following smooth step function of the gravitational permittivity was adopted: where the power index Q is the third free parameter of the theory and it regulates the transition speed between Newtonian and RG regimes (the larger its value the steeper the transition).
Being that the modification of the law of gravity is dependent on a scalar quantity, the mass density of baryonic matter 1 , it was possible to build a covariant formulation of RG [47] without the challenges encountered in defining a relativistic extension of MOND [78,[81][82][83][84]. On the other hand, a modification of the law of gravity that is dependent on a density scale might appear not so intuitive, given that the majority of the pieces of evidence on galaxy scale rather see the emergence of an acceleration scale, a 0 , below which a departure from Newtonian gravity is observed. However, the acceleration scale a 0 also seems to appear in RG, from the weak field limit (WFL) of covariant refracted gravity (CRG). This point will be better addressed in Sections 4 and 6.

Interpretation of the Three Scaling Relations in MOND and RG
The fact that the acceleration scale a 0 emerges from the three scaling relations seems to identify MOND as the most natural solution to explain them. MOND not only reproduces but actually predicted the three relations with a zero intrinsic scatter many years before they were observed, besides other pieces of evidence on the galaxy scale. It is the only theory of gravity that has this peculiar feature. Already in its first formulation of 1983 [42], Milgrom concludes that "The V 4 ∞ = a 0 GM relation should hold exactly", where V ∞ is the asymptotic flat velocity of the rotation curve, M is the total baryonic mass of the galaxy (the total mass in MOND), and the mentioned equation coincides with the BTFR (Equation (1)). The acceleration scale is set by its normalisation. Later observations confirmed this prediction (see Section 2.1). In his first works of 1983, Milgrom calculated in several independent ways the value of a 0 [42], which turned out to be consistent with the value observed years after from the scaling relations. Moreover, Milgrom predicted that the BTFR does not depend on the galaxy type or on any other galaxy property [42], in agreement with the future observations.
The condition on the acceleration a < a 0 , where the departure from Newtonian dynamics begins to be observed, can be translated in a condition on the surface mass density, Σ < Σ 0 , where Σ 0 = a 0 /G [23]. This implies a relation between the mass discrepancy in galaxies and the acceleration due to baryons, which is observed to hold from HSB to LSB galaxies with a very narrow intrinsic scatter and no dependency on galaxy properties [23].
Concerning the RAR, MOND predicted this relation with a null intrinsic scatter only for the modified inertia version of MOND and for circular orbits [79]. For the modified gravity versions of MOND, the RAR is recovered with zero intrinsic scatter only for spherical systems, whereas, for other systems morphologies, it is retrieved with a very small intrinsic scatter, in agreement with the observations [53,78,85].
A work of Eriksen et al. [86] seems to challenge MOND in modelling the RAR. Specifically, MOND might present a "cusp/core like" issue, different from the classical cusp/core problem of ΛCDM. This issue is particularly relevant for the modified inertia version of MOND. Moreover, by fitting the RAR relation (Equation (3)), in both modified gravity and modified inertia versions of MOND, from the observational data of SPARC galaxies, the best fit acceleration scale g † might result inconsistent among different galaxies.
Matsakos and Diaferio [44] demonstrated that RG reproduces both the BTFR and the MDAR. In Section 3.1.2, I showed that the asymptotic limit for the gravitational field in lowdensity environments of nonspherical systems is ∂φ/∂R ∼ a 0 |g N |, as in MOND, in the low-acceleration regime. In fact, Matsakos and Diaferio [44] write this asymptotic limit as: where b is an acceleration scale that can be expressed as: and ±h is the height from the disk plane where the condition ρ = ρ c is reached. In a simplified formulation for RG (SRG), whose conclusions can be extended to a more generic RG framework, the volume within the disk planes z = −h and z = +h is where the redirection effect of the field lines occurs.
In the low-density regime, the condition on the gravitational field given by Equation (15) can be translated in a condition on the rotation velocity: being |g N | = GM bar /R 2 . Inverting the above equation, we obtain the BTFR (Equation (1)) with the correct slope. To make the normalisation of Equation (18) in agreement with the normalisation of the observed BTFR, the acceleration b has to coincide with a 0 , which sets a condition on the z = ±h planes that depends on the galaxy baryonic mass.
Matsakos and Diaferio [44] plotted the (V f , M bar ) points from real data [87] (black dots in Figure 5) and from some disk galaxy models of star-and gas-dominated disk galaxies, with typical values of the central surface mass density, σ 0 , and of the disk scale-length, h R (open circles and squares in Figure 5). Both sets of points follow relation (18) (black solid line in Figure 5). RG also recovers the MDAR of the galaxies ( Figure 6). Matsakos and Diaferio [44] plotted the squared ratio between the rotation velocity predicted by the SRG framework and the Newtonian theory (v F and v N in Figure 6) against the Newtonian and the SRG accelerations (−N and −F in Figure 6) for the same models used to build the BTFR and for a set of point masses having m = 10 x M , with x = {7, 8, 9, 10, 11, 12}, settling between two parallel planes with z = ±h = ± √ Gm/b. The models are plotted together with a set of data points, taken from [21,23], and they properly reproduce the data, both in the (g bar , (V/V bar ) 2 ) and in the (g SRG , (V/V bar ) 2 ) planes. The gas-rich galaxies models (where the gas represents more than 20% of the total baryonic mass) have a smaller σ 0 and they distribute in the leftmost part of the MDAR, in agreement with MOND prediction according to which a larger mass discrepancy is observed when Σ < Σ 0 . On the contrary, the star-rich galaxies curves distribute on the rightmost part of the MDAR. Cesare et al. [45] built the RAR of 30 galaxies from the DiskMass Survey (DMS) [88], adopting the general RG framework and the QUMOND formulation of MOND (blue and green solid lines in Figure 7). They modelled with RG, at the same time, the rotation curve and the vertical velocity dispersion profile of each DMS galaxy, obtaining mass-to-light ratios consistent with SPS models, disk scale heights in agreement with the observations of edge-on galaxies, and RG parameters consistent among the different galaxies, suggesting their universality. To build the RG acceleration, they numerically solved the RG Poisson equation (Equation (10)), adopting the mass-to-light ratios, the disk-scale heights, and the three RG parameters estimated from the kinematic profiles of each DMS galaxy. To compute the QUMOND acceleration, they numerically solved the QUMOND Poisson equation (Equation (8)) with the same mass-to-light ratios and disk scale heights found in RG, since they are fully in agreement with the parameters found by Angus et al. [89] by performing the same analysis of the dynamics of DMS galaxies in QUMOND. The Newtonian acceleration, g bar , was calculated by solving the Newtonian Poisson equation (Equation (7)) with the same mass-to-light ratios used to compute the RG and QUMOND accelerations, being in agreement with SPS models, and with disk-scale heights h z,SR derived from the scale relation between the disk-scale lengths and heights: estimated from 60 edge-on late-type galaxies [45,90]. RG properly reproduces the asymptotic limits of the observed RAR (Equation (3), black solid line in Figure 7) but it tends to underestimate relation (3) at low g bar , even if it generally interpolates the observational data (red dots with error bars in Figure 7). Instead, QUMOND reproduces the RAR with the correct shape. This can be interpreted by the fact that RG might attribute more luminous mass than QUMOND.
A more serious problem for RG is that the RAR presents a too large intrinsic scatter (0.11 dex, whereas the possible intrinsic scatter of the RAR found by Li et al. [56] for SPARC galaxies is 0.057 dex) and some correlations between the residuals from Equation (3) and certain galaxy properties, in disagreement with the observations [22,53]. However, this question has to be further deepened by building the RAR in RG for a larger sample of disk galaxies with more accurate rotation curves, such as SPARC, before concluding that this is due to an issue of RG theory. Indeed, the RAR of DMS data also shows some correlations between the residuals and some galaxy properties, which might not be observed in the SPARC sample and suggests that the DMS could not be the most suitable sample where to investigate the RAR.
In contrast, the RAR curves computed in QUMOND, very neatly distribute around Equation (3) with an intrinsic scatter of 0.017 dex. This is consistent with expectations, since QUMOND is a modified gravity version of MOND, and, thus, it does not provide a RAR with zero intrinsic scatter for nonspherical systems as disk galaxies [78]. Moreover, QUMOND RAR also presents correlations between its residuals and some galaxy properties, again in agreement with the fact that, in this case, the scatter of the RAR is not equal to zero and which might further suggest that the DMS sample was not the most suitable to investigate the RAR.

Possible Interpretations for an Intriguing Acceleration Scale
If the observed DM-baryons scaling relations might appear not so intuitive in the ΛCDM context, a piece of evidence even more difficult to interpret in this framework is the emergence of an acceleration scale, a 0 , from the three relations. However, this acceleration scale seems to be retrieved in the ΛCDM model as well, for example, from the MDAR and the RAR built from the simulated galaxies in the Magneticum simulation at z ∼ 0.1 [54]. These relations were fitted with Equations (2) and (3), obtaining an a 0 consistent with 1.2 × 10 −10 m s −2 [54] (see Section 2.2). Mayer et al. [54] also found that Magneticum galaxies followed a MDAR (Equation (2)) and a RAR (Equation (3)) relation at higher redshifts but with a fitted a 0 substantially different. Specifically, a 0 decreases for decreasing redshift, that is, in the more recent Universe. This trend means that the mass discrepancy decreases, and, thus, galaxies become more baryonsdominated, as cosmic time advances, and this can be explained by the progressive cooling of the gas during time.
Mayer et al. [54] tried to assess whether the a 0 (z) relation observed in simulated Magneticum galaxies is consistent or not with MOND predictions, which is not trivial since MOND is not formulated as a relativistic theory. Milgrom [41] pointed out that a 0 ≈ cH 0 , which might suggest that the redshift dependence of a 0 is similar to the one of the Hubble parameter H(z). This would imply the following relation: where the fraction of baryonic matter Ω m and of the cosmological constant Ω Λ with respect to the total mass-energy budget of the Universe are put equal to 0.25 and 0.75, following [91], and a 0 (0) is the value of the acceleration scale a 0 = 1.2 × 10 −10 m s −2 observed today. However, this relation is too steep compared to the a 0 (z) relation observed in Magneticum simulation and it also seems to disagree with observational data. More realistic predictions of the a 0 (z) relation in MOND might be provided by relativistic theories that have MOND as a limiting case in a nonrelativistic regime, such as tensor-vector-scalar gravity (TeVeS) [92] and covariant emergent gravity (CEG) [93]. TeVeS might suggest that a 0 varies on timescales larger than the Hubble time, even if, according to [91], this change is also possible on cosmological timescales. Anyway, the a 0 (z) relation predicted by TeVeS is not consistent with the one observed in [54] from Magneticum galaxies. CEG predicts an a 0 (z) relation dependent on the variation of the size of the cosmological horizon [94]. However, the redshift dependence of this a 0 (z) relation is smaller than the one observed from the galaxies in the MUGS2 simulated sample in ΛCDM [71]. An even more striking coincidence is given by the near equality between a 0 emerging from the local universe, that is, from galactic dynamics, and from cosmology. Given the two acceleration cosmological parameters [95]: and where H 0 and Λ are the present values of the Hubble parameter and of the cosmological constant, the following equivalence subsists: This arises an outstanding coincidence that needs careful investigation since Λ is assumed to be constant across cosmic time and H 0 is the value at the present epoch of the Hubble parameter, H(t), that varies across cosmic time.
Some interpretations of the relation a 0 ∼ Λ 1/2 in the MOND context were provided by [23,[95][96][97]. In the most recent work among these, Milgrom [95] re-wrote Equation (23) in terms of length or mass, introducing the "MOND length", l M , and the "MOND mass", M M : and where l H ≡ cH −1 0 is the Hubble radius, l Λ ≡ (Λ/3) −1/2 is the de Sitter radius related to Λ, and M H and M Λ are the total mass of the Universe enclosed within l H and l Λ , respectively.
Equation (23) emerges in some particular effective-field, relativistic extensions of MOND, where MOND is retrieved in their WFL. In these MOND relativistic formulations, an additional Lagrangian term, is added to the GR Einstein-Hilbert (EH) Lagrangian density: where R is the Riemann curvature scalar. In Equation (26), l is a length constant, needed to provide the Lagrangian with the correct dimension, F is a dimensionless function, and Q has the dimension of a length −2 and it is built from the first space-time derivatives of the gravitational degrees of freedom [95]. Some examples of MOND relativistic extensions where a Lagrangian-like term as Equation (26) is included are the MOND adaptations of the Einstein-Aether theories [98], bimetric MOND (BIMOND) [99], bimetric massive gravity [100], and the noncovariant theory presented in [101]. If we include the cosmological constant term in the EH Lagrangian, Equation (27) becomes: Therefore, we can see from Equation (26) that any constant term added to F (x), which has to be of the order of unity due to naturalness, can be identified with a cosmological constant term of this kind: that, combined with Equation (24), gives the a 0 -cosmology coincidence given by Equations (22) and (23) [95]. It is important to point out that this result is not obtained by adding a cosmological constant term, Λ, "ad-hoc" to the Lagrangian, verifying a posteriori that Λ ∼ (a 0 /c 2 ) 2 , but by the presence of the two l-terms in the Lagrangian (26), appearing inside and outside F (x) for dimensional reasons and playing different roles. This suggests that the two l-terms derive from the same underlying physics, whereas a Λ constant added "ad-hoc" might have suggested a different physical origin. Milgrom [95] also suggested that the a 0 -cosmology connection can also emerge from a scenario in which the Universe is seen as a sphere-like submanifold, that is a brane, of radius l H or l Λ , embedded in a space-time with higher dimension [102]. The dynamics of the submanifold is the one we observe, that is the MOND one, and it emerges from the dynamics of the higher dimension space-time. In fact, the acceleration scale a 0 is seen as an emergent acceleration constant in a brane scenario. Further details can be found in [102].
In RG, where the transition between Newtonian and modified gravity regimes is regulated by a density scale ρ c , rather than by an acceleration scale a 0 , the a 0 -cosmology relation might appear harder to interpret. However, Sanna et al. [47] recently formulated a covariant extension of RG, where a 0 seems to emerge from the WFL of the theory and the a 0 -Λ coincidence might be explained.
CRG is a scalar-tensor theory that introduces a single scalar field ϕ, nonminimally coupled to the metric, that mimics both the DM effect on galaxy scale and the DE effect on cosmic scale, namely, the accelerated expansion of the Universe. Therefore, CRG belongs to the restricted class of modified gravity models that invoke a unified dark sector, that is, that attribute the phenomenologies of DM and DE to a single cause.
CRG is derived from a general scalar-tensor action: (see [47] for the explanation of the symbols). In CRG, the general differentiable function of the scalar field W (ϕ) is: and the potential V (ϕ) has a self-interaction form: where Ξ is a constant parameter. With the definitions adopted for W (ϕ) and V (ϕ), the CRG equations below are obtained: and CRG is based on a chameleon screening mechanism, e.g., [103], that is, in regions where the Newtonian WFL holds, the extra degree of freedom of the scalar field ϕ mediates a fifth force which can be detected, whereas, in high-density regions, this degree of freedom is screened. This behaviour is also what we expect from the RG gravitational permittivity (ρ), which, therefore, might be related to ϕ.
The WFL of CRG holds the following equation: which reduces to RG Poisson equation (10) if the scalar field is twice the permittivity, ϕ = 2 . This is an important result since it confirms that the scalar field is associated with the phenomenology mimicked by the gravitational permittivity, that is of DM on galaxy scale. Instead, the Newtonian Poisson equation is recovered for a constant scalar field ϕ = 2.
Calculating the CRG gravitational field in the WFL from Equation (35) for a spherical source with density ρ s (r), monotonically decreasing with r, immersed in a homogeneous background with constant density ρ bg , Sanna et al. [47] found that an acceleration scale can be set. At large distances from the source, the scalar and the gravitational fields, ϕ and dφ/dr, are linked by the relation: where ρ(r) is ρ s (r) + ρ bg . The acceleration scale: is set from Equation (36). In regions where dφ/dr a Ξ , the gravitational field has a similar rdependence as the gravitational field calculated close to the source, that is, the Newtonian field. Instead, for dφ/dr a Ξ , it departs from the Newtonian one. From this result, the acceleration a Ξ recalls a 0 since it demarcates Newtonian from modified gravity regimes. Solving the CRG field equations for a homogeneous and isotropic universe with flat curvature, Sanna et al. [47] found that Ξ ∼ Λ, and, thus, Ξ plays the role of the cosmological constant in ΛCDM. In the limit obtained at large distances from the source, Equation (37) becomes: Using the observed value of Λ at the present epoch, Equation (39) holds a Ξ ∼ 10 −10 m s −2 , fully in agreement with the value of a 0 , and it also provides the relation a 0 ∼ Λ 1/2 , which is the observed a 0 -cosmology connection.
The difference between a Ξ and a 0 consists in the fact that, whereas a 0 is a constant independent of the gravitational field source, ρ, a Ξ depends on the source, given Equation (37), even if for large distances from the source this dependence drops since ρ s ρ bg . Future investigations have to verify whether by repeating the calculations for a generic case and not for a specific source, the connection between a Ξ and a 0 continues to hold in a real Universe.
Given the connection between Ξ and the cosmological constant Λ, it can be concluded that the scalar field ϕ is also responsible for the accelerated expansion of the Universe, that is for the DE phenomenology, besides the DM phenomenology on galaxy scale, given its connection to the gravitational permittivity . The fact that the DM and the DE effects are described by a single scalar field might provide an advantage for CRG since the idea of a unified dark sector is theoretically justified. The a 0 ∼ Λ 1/2 relation itself represents one possible piece of evidence for the unification of the dark sector. Moreover, Martin Kunz claims that gravity can only probe the total energy-momentum tensor of the Universe, implying a degeneracy between the dark components, which further goes in this direction [104,105].
Among these models, Conformal gravity is one of the first to invoke the unified dark sector, even if it might present some problems in modelling the rotation curves of galaxies and the gravitational lensing effects [108]. Two recent and novel models invoking the unified dark sector are the unified superfluid dark sector [111] and the fuzzy dark fluid [125]. In the unified superfluid dark sector, the Universe is dominated by a unique DM superfluid made of axion-like particles, with two energy states having an energy gap smaller than H 0 that can interact with each other. These interactions at the microscopic level change the macroscopic behaviour of the fluid, producing an accelerated expansion of the Universe that mimics DE. Besides the effect of DE, this fluid can also mimic the galaxy phenomenology due to DM without facing some problems encountered by others of these models, such as superluminal sound speeds or the need for a UV completion.
In the fuzzy dark fluid model, besides the behaviours of DM and DE, a single scalar field also mimics the behaviour of inflation, by assuming a nonminimal coupling to the gravitational field, a Mexican hat-shape potential, and a spontaneous symmetry breaking before the inflationary period. This peculiar feature of a unique description of the DM, DE, and inflation phenomenologies is also shared by mimetic gravity [123].

Dwarf Galaxies and Globular Clusters
The flatness of the rotation curves of disk galaxies and the three scaling relations, BTFR, MDAR, and RAR, observed on the galaxy scale, are only some of the predictions of MOND. MOND predicted additional pieces of evidence about galactic dynamics before they were observed. Among them, the dynamics of LSB galaxies is worth mentioning.
As already mentioned in the Introduction (Section 1), HSB disk galaxies are dominated by stellar mass in their central regions, where their rotation curves are steeply rising toward their asymptotic values. Instead, the dynamics of LSB galaxies, generally dwarf and dSph galaxies, is observed to be different. Their rotation curves are slowly rising toward their flat region and they appear to be DM-dominated even in their innermost regions. Therefore, the maximum-disk hypothesis cannot be applied as for HSB galaxies [31,32]. Specifically, dwarf galaxies are among the known darkest galaxies observed: they have an inner velocity dispersion σ ∼ 10 km s −1 , an order of magnitude larger than the velocity dispersion σ ∼ 1 km s −1 expected for systems having the same luminosity and scale radius (∼100 pc) at equilibrium [33]. Their luminosity varies in the range ∼[10 2 , 10 10 ] L [126,127] but their velocity dispersions are similar, which might indicate that they are dominated by a similar DM distribution [19,128].
This different shape of the rotation curves of HSB and LSB galaxies reflects in one of the small-scale issues of the ΛCDM model, the cusp/core problem. Indeed, steeply-rising rotation curves can be modelled by a Navarro Frenk White (NFW) DM density profile [129,130], cuspy in its innermost part and predicted by collisionless N-body simulations, whereas slowlyrising rotation curves can only be reproduced by a cored DM density profile, which might be accounted for in ΛCDM only introducing baryonic feedback and tidal effects [131].
The different dynamic properties of these two categories of galaxies were instead predicted by MOND some years before they were observed [42]. Milgrom [42] predicted that dwarf galaxies would have shown strong deviations from standard gravity and in particular that when the velocity dispersion data of dwarf galaxies had been available, these galaxies would have presented a mass discrepancy equal to 10 or larger, depending on their distance from the Milky Way, when modelled in standard gravity.
As already anticipated in Section 3.2, the MOND acceleration scale a 0 can be translated on a surface mass density scale Σ 0 = a 0 /G [23,42]. Small surface mass densities also indicate small surface brightnesses and Milgrom predicted that LSB galaxies would have shown stronger effects of the modification of the law of gravity with respect to the Newtonian one. The effect of gravity modification already appears in the innermost radii of these galaxies. Milgrom [42] predicted that a transition radius r 0 between the standard and the modified gravity regimes, dependent on the local value of the rotation velocity V, would have been set when the equality V 2 /r 0 ≈ a 0 occurs. In particular, where V 2 /r a 0 , the local mass-to-light ratio of the galaxy should not indicate the presence of a hidden mass, and where V 2 /r starts to go below a 0 , the local M/L should begin to rapidly increase. The smaller the average surface brightness of the galaxy, the smaller the r 0 in units of the galaxy scale length h R . A similar result was also found by [44] for RG (see the top panel of Figure 10 of [44]). Moreover, Milgrom [42] also predicted a correlation between the average surface mass density or surface brightness and the steepness of the rotation curve in reaching its asymptotic limit, which states that galaxies with small surface densities show a slowly-rising rotation curve and vice-versa, as observed by successive measurements.
A possible challenge for MOND on galaxy scale is instead provided by the internal dynamics of GCs, which have baryonic masses similar to the ones of LSB galaxies, settling in the outermost regions of the Milky Way. In those regions, the background acceleration is much smaller than a 0 and the external field effect is negligible. Therefore, we expect that their stellar velocity dispersion profiles present a MONDian behaviour but this is not the case. Newtonian theory without the presence of a DM halo can better fit these kinematic profiles than MOND, whose predicted velocity dispersions can exceed the Newtonian ones also by a factor of ∼3 [34][35][36][37][38][39][40]. Relevant examples for this result are provided by the GCs NGC 2419, Palomar 14, and Palomar 4 [35,36,[38][39][40]]. Yet, this tension between MOND predictions and measurements might not indicate an issue of the theory but it might be due to inaccurate data or approximate modelling.
Concerning the former case, inaccurate data can derive from low-resolution spectroscopy and from errors on GCs distances larger than 10% of their values (e.g., [37]). The Gaia mission, which provides the parallaxes of the stars from which to derive their distances, accurate at the µarcsec level, might represent a turning point in this sense [132]. Concerning the latter case, most of the adopted models assume spherical symmetry, absence of rotation, and orbital isotropy. Whereas the first two assumptions can be justified by observations, orbital anisotropies in GCs are predicted by N-body simulations [133] and including strong radial anisotropies in the modelling can reconcile MOND expectations with the observed velocity dispersions. Further studies with high-precision measurements, which allow us to neatly disentangle the effect of a strong radial anisotropy and the adopted theory of gravity, need to be performed [37]. Moreover, the question might be even harder to understand, since some theories of formation and evolution of GCs predict the presence of DM in these objects and its observational evidence is debated [134][135][136].
Instead, RG theory provides a more natural solution for the different dynamic properties of LSB, dwarf, and dSph galaxies and GCs. As anticipated in Section 3.1.2, RG predicts a different shape for the gravitational field lines in flat and spherical systems. Specifically, whereas in spherical systems the gravitational field lines always maintain a radial direction following the Newtonian trend and become enhanced compared to the Newtonian field in the external regions of the systems, where the density goes below the critical value ρ c , in nonspherical configurations we observe a refraction of the field lines toward the equatorial plane of the object and the deviations from Newtonian gravity already begin to appear in regions where ρ > ρ c . In particular, the flatter the system the stronger the redirection effect of the field lines and the larger the mass discrepancy if interpreted in a Newtonian framework. This can intuitively explain the diverse dynamic behaviour of LSB, dwarf, and dSph galaxies and GCs, the former generally having a flatter shape and the latter a more spherical one [44].
This RG prediction is also in agreement with a claimed positive correlation between the elliptical galaxies' ellipticities and their DM content [137,138].
As mentioned in Sections 1 and 3.2, Cesare et al. [45,46] demonstrated that RG can model the dynamics of flat (30 disk galaxies from the DMS) and spherical (three E0 galaxies from the SLUGGS survey) systems. The modelling of the two classes of systems is obtained with statistically consistent { 0 , Q, ρ c } parameters, showing that the gravitational permittivity is independent of the shape of the considered system. To perform a more robust test of the theory, Cesare et al. [45] estimated the three RG parameters both from each individual DMS galaxy, simultaneously modelling its rotation curve and vertical velocity dispersion profile, and from the kinematic profiles of the entire DMS sample considered at the same time. The two sets of values are consistent within 2σ. The average Q and ρ c derived from the simultaneous modelling of the velocity dispersions of the stars and the blue and red GCs in each E0 galaxy are consistent within 1σ with the Q and ρ c averaged from the values obtained from each DMS galaxy, whereas the 0 parameters are still in agreement within 3σ. The average Q and ρ c from the E0 galaxies are also consistent, within 3σ, with the unique combination of Q and ρ c derived from the entire DMS sample. Instead, the tension increases to 14.8σ for the 0 parameter. However, this does not necessarily indicate an issue for RG. This might be due to the approximate procedure with which the single combination of RG parameters is estimated from the entire DMS sample, which results in error bars of the 0 parameter that are much smaller than the error bars of the average 0 of DMS galaxies. Moreover, it might be due to incorrect modelling of the dynamics of elliptical galaxies, which are treated as isolated systems without net rotation, or to a wrong assumption for the gravitational permittivity functional form [46]. A review of both the works about disk and elliptical galaxies in RG is presented in [139].

Discussion and Conclusions
The most investigated cosmological model is ΛCDM, which assumes the validity of GR and the inclusion of two dark components, DE and DM, besides baryonic matter, which only represents the ∼5% of the mass-energy budget of the Universe. The ΛCDM paradigm reconciles with the majority of the observations, from the largest to the smallest scales. However, the results of the detection, through direct, indirect, or collider experiments, of the most investigated DM candidate, the weakly interacting massive particles (WIMPs), are still under debate [140]. Moreover, the nature of DE is still unknown. Future experiments, such as Euclid, might shed light on its nature.
Furthermore, the ΛCDM model presents some issues, both on large and on small scales. Particularly relevant are the problems observed on the galaxy scale, such as some remarkable coincidences that can hardly be explained by the stochastic merging process of structure formation predicted by ΛCDM, unless precise fine-tuning between DM and baryonic processes is invoked. Among these coincidences, three scaling relations between dark and baryonic matters in galaxies, the BTFR, the MDAR, and the RAR, that neatly quantify the mass discrepancy on galaxy scale, are observed and they see the appearance of the acceleration scale a 0 = 1.2 × 10 −10 m s −2 . The results of Mayer et al. [54] from the galaxies in the Magneticum simulation can reproduce in ΛCDM the emergence of this acceleration scale at redshift z ∼ 0.1. They also predicted an evolution of a 0 with redshift, where a 0 decreases with increasing redshift, which still has to be confirmed. Intriguingly, the acceleration scale a 0 presents another coincidence, being its value consistent with the combination of some cosmological parameters: a 0 ∼ H 0 ∼ Λ 1/2 . In particular, the relation a 0 ∼ Λ 1/2 links a quantity observed on galaxy scale and the parameter that regulates the Universe accelerated expansion, which suggests a unification of the two dark sectors and connects the physics on small and large scales. This is even less intuitive to interpret in a ΛCDM framework.
On the other hand, MOND theory not only reproduces but even predicted, with either a small or a null intrinsic scatter, these three relations, assuming a modification of the law of gravity for accelerations smaller than a 0 . The value of a 0 was estimated by Milgrom [42] in several independent ways before it emerged from observations and it turned out to be consistent with the value of a 0 observed some years later from the DM-baryons scaling relations. MOND predicted other pieces of evidence on the galaxy scale, such as the fact that LSB galaxies appear "darker" and with more slowly-rising rotation curves than HSB galaxies. However, it presents some issues in describing the dynamics of GCs residing in the Milky Way outskirts, which present a Newtonian behaviour even if the background acceleration is below a 0 . MOND also reproduces the a 0 ∼ Λ 1/2 relation, as presented in several studies [23,[95][96][97].
A more recent theory of modified gravity, RG, formulated in a nonrelativistic way by Matsakos and Diaferio in 2016 [44], might shed further light on galaxy dynamics. RG has already presented some encouraging results on galaxy scale, reproducing the dynamics of both flat and spherical systems with a consistent set of RG parameters [45,46,139], the BTFR, and the MDAR of real and simulated galaxies. It also models the RAR of DMS galaxies, even if with some issues, which requires further investigations to assess whether these problems are due to the chosen galaxy sample or to RG theory. RG predicts a different shape for the gravitational field lines in spherical and nonspherical systems. Specifically, the field lines remain radial in spherical systems and become increasingly refracted toward the equatorial plane of increasingly flat systems. The refraction of the field lines produces a boost of the gravitational field that mimics the presence of a DM halo in Newtonian gravity. This means that an increasingly flat system is increasingly DM-dominated, if interpreted in the Newtonian context. With this feature, RG naturally due to the different dynamic properties of LSB galaxies and GCs to their diverse shape, the former generally being flatter and the latter nearly spherical. This system morphology-mass discrepancy relation predicted by RG is also consistent with the possible ellipticity-total M/L correlation of elliptical galaxies estimated by Deur [137,138].
Despite the promising results shown by different theories, we are still far from the answer for a final scenario of the cosmological model. Given the emergence of a 0 from several pieces of evidence on galaxy scale, MOND might seem the most intuitive solution. However, it presents several issues on a larger scale. For example, it can reduce but not eliminate all the mass discrepancy in galaxy clusters [141][142][143]. Moreover, the building of a covariant version of MOND seems to appear challenging. Some attempts at formulating this relativistic extension failed to describe the features of gravitational lenses, provided superluminal speeds, or were not in agreement with the post-Newtonian tests of General Relativity [78,81,82]. A relativistic extension of MOND is TeVeS [92], which solved some of these problems but was unable to reproduce cosmological pieces of evidence, such as the CMB or the matter power spectra [83,84]. However, further studies about covariant MOND are still ongoing and some recent results might look promising [144,145].
RG might provide an alternative solution. Despite being based on a density scaledependent modification of the law of gravity, which is not what observations might suggest, a covariant formulation of RG, CRG, seems to be promising, given the results of [47] that also show that the acceleration scale a 0 might emerge from the WFL of the theory. In particular, building a relativistic extension of RG was possible since the modification of the law of gravity depends on a scalar quantity, in this case the density, whereas for MOND it depends on a vector quantity, namely, the acceleration. CRG describes the DM and DE phenomenologies with a single scalar field, suggesting a unified dark sector, and retrieves the a 0 ∼ Λ 1/2 relation, which is a remarkable result.
However, further studies have to be performed to validate RG. A more accurate study of elliptical galaxies, removing the assumptions adopted in [46] and considering a larger sample with different ellipticities and extended kinematic profiles (e.g., SLUGGS [146] and ePN.S [147] surveys), have to be made to better assess if RG can reproduce the dynamics of these systems. The fact that RG can account for the different dynamics of dwarf galaxies and GCs is only a hypothesis that should be tested on real samples, such as the dwarf galaxies surrounding the Milky Way, e.g., [148] and belonging to the LITTLE THINGS survey [149], and the GCs settling in the Milky Way outskirts, in particular NGC 2419, Palomar 14, and Palomar 4, e.g., [35,36,39,40]. Moreover, RG should be tested on larger scales, to verify whether it can describe the dynamics of galaxy clusters. Some preliminary encouraging results in this sense were obtained by Matsakos and Diaferio [44] but these studies have to be extended to larger samples (e.g., CIRS and HeCS [150,151]) and with less approximate modelling. At last, the studies on cosmological scales have to be completed with CRG.
MOND and RG are only two possible theories of modified gravity but many other theories that might provide a solution for small-and large-scale problems have been built. Another theory that is worth mentioning is scalar-vector-tensor gravity (SVTG), better known as modified gravity (MOG) [152]. Whereas the modification of the law of gravity in MOND and RG depends on an acceleration and on a density scale, respectively, in MOG it depends on a length scale. MOG is a theory of gravity built in a covariant way that introduces a scalar, a tensor, and a massive vector field, whose contributions are added to the classical EH action. MOG assumes that the gravitational constant G, the coupling constant ω and the mass µ of the vector field, and the cosmological constant Λ are dynamical scalar fields which vary with space and time [152].
MOG has two progenitor theories, nonsymmetric gravity theory (NGT) [153] and metricskew-tensor gravity (MSTG) theory [154], which produce the same modified acceleration law as MOG for weak gravitational fields. NGT, MSTG, and MOG presented several encouraging results on different scales. They can reproduce solar system and terrestrial gravitational tests, the observations of the binary pulsar PSR 1913+16 [155], the rotation curves of both HSB and LSB galaxies and the BTFR [154][155][156], the dynamics of an elliptical galaxy [155], the velocity dispersion of Milky Way (MW) satellite galaxies [157], the internal velocity dispersion profiles of GCs [158], the cluster lensing [154], the mass profiles of galaxy clusters derived from X-ray emitting gas [159] and their thermal profiles [156], the Bullet cluster [160], and several cosmological observations, such as the CMB temperature anisotropy, the galaxy power spectrum, and the supernova luminosity-distance measurements [156,[161][162][163]. A more recent work seems to demonstrate that MOG is able to reproduce also the RAR of 153 spiral galaxies from SPARC sample [164].
The dynamics of GCs, of galaxy clusters, and of systems at sufficiently large distances from their centres provide an important test to distinguish among MOG, MOND, and RG. MOG reproduces the internal velocity dispersions of GCs around the MW, independently from their distance from the Galaxy centre, consistently with Newtonian expectations. This is not the case for MOND, which predicts a MONDian behaviour for the dynamics of GCs sufficiently distant from the MW centre, such that the background acceleration goes below a 0 . However, observations suggest a Newtonian behaviour also for these GCs, where MOND predictions may exceed the measured velocity dispersions also by a factor of ∼3, e.g., [35]. Moreover, as already specified, whereas MOG can reproduce the masses of galaxy clusters, and also RG is likely to do this, MOND can only retrieve them with an additional DM component.
To distinguish among MOG, MOND, and CRG, it would be ideal to have the kinematic data of different systems at large distances from their centres. Indeed, at these distances, MOG predicts a Newtonian keplerian trend of the rotation velocity, i.e., V(R) ∝ R −1/2 , consistent with the results of [165] and with the gravitational lensing results of [166]. However, this trend is not predicted by MOND and RG (see Section 3).
Both MOG and RG, also with its covariant extension CRG, reproduce pieces of evidence on different scales. However, some fundamental features distinguish MOG from CRG. Whereas MOG introduces a scalar, a tensor, and a vector fields, CRG only includes a single scalar field. Moreover, whereas the CRG scalar field is responsible for both the DM and the DE phenomenologies, MOG attributes the two dark sectors to two different causes and it introduces a cosmological constant, dependent on space and time, to mimic DE effect. Concerning the a 0 ∼ Λ 1/2 relation, whereas CRG retrieves it as a consequence, MOG imposes it [155] to constrain some MOG parameters (G 0 , M 0 , and r 0 ). Several works, e.g., [155,156], demonstrate that MOG can model different pieces of evidence with a minimal number or with no free parameters. The results of Cesare et al. [45,46] show that a universal combination of RG free parameters might exist. Additional studies have to be performed to verify if the phenomenology at several scales in the Universe can be modelled with this unique combination of RG free parameters. Generally, further tests have to assess which of the three theories of modified gravity, MOND, RG, or MOG, better describe the pieces of evidence on different scales.
To conclude, the correct cosmological model might, at last, be ΛCDM and further baryonic processes might still be discovered to properly reconcile the theory with all the observed coincidences on the galaxy scale.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: In fact, the gravitational sources depend on other scalar quantities besides the mass density, such as their total mechanical and thermodynamical energy or their entropy. Yet, these quantities depend in turn on the mass density and, thus, adopting a gravitational permittivity also dependent on these quantities would be likely to produce a phenomenology comparable to the one obtained with the simple dependence on the mass density alone [44,45].