Modified Newtonian Gravity, Wide Binaries and the Tully-Fisher Relation

Abstract

A recent study of a sample of wide binary star systems from the Hipparcos and Gaia catalogues has found clear evidence of a gravitational anomaly of the same kind as that appearing in galaxies and galactic clusters. Instead of a relative orbital velocity decaying as the square root of the separation, $\Delta V\propto r^{-1/2}$, it was shown that an asymptotic constant velocity is reached for distances of order $0.1$ pc. If confirmed, it would be difficult to accommodate this breakdown of Kepler’s laws within the current dark matter (DM) paradigm because DM does not aggregate in small scales, so there would be very little DM in a $0.1$ pc sphere. In this paper, we propose a simple non-Newtonian model of gravity that could explain both the wide binaries anomaly and the anomalous rotation curves of galaxies as codified by the Tully-Fisher relation. The required extra potential can be understood as a Klein-Gordon field with a position-dependent mass parameter. The extra forces behave as $1/r$ on parsec scales and r on Solar system scales. We show that retrograde anomalous perihelion precessions are predicted for the planets. This could be tested by precision ephemerides in the near future.

1. Introduction

In the 1930s, the astronomer Fritz Zwicky conducted pioneering work in which he examined the validity of the virial theorem for the Coma galaxy cluster, unveiling a surprising anomaly. Zwicky discovered that galaxies within the cluster move much faster than predicted by the gravitational attraction of the observed luminous matter [1]. Famous for his bold proposals, he interpreted this observation as the consequence of a vast amount of matter that cannot be directly observed in the cluster. This matter was called “dunkle (kalt) materie” or dark matter. Although Zwicky’s idea was dismissed and ignored by his contemporaries, there is no doubt that this early interpretation has influenced modern thoughts on this—still unresolved—problem in astrophysics and cosmology. In 1936, Smith replicated Zwicky’s work to estimate the mass of the Virgo cluster. He also concluded that a great mass of intergalactic material should be present within the cluster [2].

Rubin et al. [3] and Faber & Gallagher [4] found similar results for individual galaxies in the late 1970s and early 1980s. They showed that the rotation curves of spiral galaxies reach a “plateau” at large distances from the center instead of decreasing roughly as $1/r^{1/2}$, where r is the distance from the center as implied by Kepler’s laws.

Later on, the concept of dark matter has been incorporated in the—now standard—$\Lambda$CDM model of cosmology [5,6,7,8,9]. This concept has found a place in the cosmological models because of its explanatory power but also because there has been no widespread interest in alternative explanations.

In the absence of DM, cosmological models would predict a universe with an age at odds with that inferred from the oldest objects it contains. Metal-poor stars are among the first stars formed after the Big Bang, consequently they are characterized by a low concentration of iron and other heavy elements. These elements were formed in core-collapse supernova explosions whose initial proportion can be estimated. From the relative abundance of Thorium and Uranium isotopes in these metal-poor stars, very good lower bounds on the age of the Universe have been obtained [10,11,12]. An age of $14\pm 2$ Gyr was deduced from these studies, concordant with the age deduced from the $\Lambda$CDM model with a fraction ${\Omega }_{DM}=0.23$ of dark matter. The age of globular clusters can also be inferred from the turning point in their Hertzsprung-Russell diagrams. This is also in agreement with the age of the Universe in the standard $\Lambda$CDM model [13]. Without DM, the primordial abundance of ${}^4$He would also be much larger than the accepted value of about 25%. If all the missing matter would be baryonic, nucleosynthesis would have begun earlier (because of the higher reaction rates) and, consequently, at an era with a higher fraction of still undecayed neutrons. These should have ended up as deuterons and later as ${}^4$He [14]. As this scenario is not strongly disfavoured observationally, we must conclude that DM, if it exists, is mainly, non-baryonic.

The most popular candidate for the DM particles are WIMPs (Weakly Interacting Massive Particles) with masses in the range of 10–${10}^3$ GeV/c² [15]. Tight constraints have been, experimentally, imposed on the cross-sections for their weak interactions with nucleons [16,17]. In the DAMA/LIBRA experiment, an annual modulation was detected that was considered to arise from the interaction of DM particles in the Galactic halo with the detector mass as the velocity of the Earth varies along its orbit. However, this modulation has been fitted using conventional interactions with muons and neutrinos [18]. Moreover, other experiments have not been able to reproduce the results of DAMA/LIBRA after adequate corrections for surface effects and the expected background [19,20].

Apart from the lack of solid evidence about the existence of WIMPs from direct detection in underground experiments, there are other more general objections against the hypothesis of cold dark matter (CDM). Structure formation models using CDM are characterized by two main problems as comparisons with observations of galaxies and clusters are carried out: (i) Most dwarf galaxies show a flat density profile of DM at their core (as inferred from their rotation curves). On the other hand, simulations predict that a DM density increasing at small radii (in the core). This discrepancy among theory and observations is known as the “cuspy halo problem” [21]; (ii) Hierarchical clustering of dark matter during the formation of clusters yields a large number of dwarf satellite galaxies as predicted by simulations. The real number of satellite galaxies is too small, at odds with the simulation results [22].

These failures of the standard CDM model have been strongly criticized by Kroupa and collaborators [23,24,25]. These authors have made a convincing thesis concerning the necessity of a new paradigm in cosmology beyond the dark matter hypothesis. Another important issue is the so-called planes of satellite galaxies problem as reviewed by Pawlowski [26]. It has been shown that the satellite galaxies of the Milky Way and of Andromeda galaxy align in almost flattened planes but these anisotropic distributions are extremely rare in $\Lambda$CDM simulations of galaxy formation. In the case of the Andromeda galaxy, the proper motion measurements also suggest that the satellite members are co-rotating within a plane [27]. By proposing a flyby of the Milky Way and the Andromeda galaxy $\simeq 8$ Gyr ago this particular distribution can be explained in the context of MOND as shown by Banik et al. [28]. This issue is very important because the only known explanation is that the Local Group and Centaurus A satellite planes are tidal dwarf galaxies, which should be free of DM in the CDM picture. The problem with this idea is that these satellites have high internal velocity dispersions impliying a breakdown of Newtonian dynamics because we cannot use the standard explanation in terms of CDM.

Other proposals for a dark matter particle candidate such as axions [29,30], magnetic monopoles [31,32] and supersymmetric partners of known particles [33,34] share similar difficulties.

Another candidate that has attracted much attention in recent times is the sterile neutrino [35,36]. These sterile neutrinos may have formed by oscillations in the early Universe [37], lepton-number-driven resonant conversion [38] or decay of a heavy scalar [39,40,41,42]. It is supposed that they have a mass of a few keV/c² making them a form of “warm dark matter” (WDM). This term implies that the particles presently move at velocities similar to that of stars and gas clouds within galaxies. On the other hand, their spectra can be very different from a Maxwellian distribution [36]. Sterile neutrinos might be a way to overcome the cuspy-halo problem [35].

Evidence for the existence of a sterile neutrino sector beyond theoretical speculations is, however, small. The discovery of a $3.5$ keV line in the X-ray spectra of some galaxies has been interpreted as originating from the decay of a 7 keV/c² sterile neutrino. This interpretation is—surely—too hasty and premature. Gu et al. [43] suggested a conventional mechanism in atomic physics (charge exchange with bare sulfur ions) that can explain the generation of this particular line. Other authors have found no evidence of such a line in stacked galaxy spectra [44] or from our own Galaxy, which should be emitting such X-rays from all directions of the Galactic halo [45].

The lack of experimental evidence for DM particles and the inconsistencies of the DM paradigm have motivated throughout the years several alternative models and theories. These models share the common idea that some modification is required to Newtonian dynamics, Newtonian gravity and General Relativity. These modified gravity theories are expected to explain away DM as a “phantom mass” not corresponding to any real particle but to a breakdown of our current gravity theories on galactic and intergalactic scales, and perhaps also in the Solar system. One of the earliest models of this kind is M. Milgrom’s modified Newtonian dynamics (MOND) [46,47,48,49]. MOND hypotheses lead to a modification of Newtonian gravity:

$$g_N=g\mu \left({\displaystyle \frac{g}{a_0}}\right),$$

(1)

where $g_N$ is the Newtonian gravity and g is the true gravity. Here $\mu \left(x\right)$ is a function with the limits $\mu \left(x\right)\to 1$ as $x\to \infty$ and $\mu \left(x\right)\to x$ as $x\to 0$. With the parameter $a_0$ tuned to the value $a_0=1.2\times {10}^{-10}$ m/s${}^2$, MOND can explain both the anomalous rotation curves of galaxies and the Tully-Fisher relation among luminosity/baryonic mass of galaxies and asymptotic rotation velocity [49,50]. Moreover, MOND provides a solution to the “cuspy-halo problem” because for the core region of galaxies $a>>a_0$ and the standard Newtonian behaviour is recovered. MOND also has its own problems, the most conspicuous being the need for some extra mass in galactic clusters to obtain agreement with observations [51]. One possible solution to this cluster scale problems of MOND is a hybrid MOND+hot dark matter model as proposed by Angus [52]. It has also recently been shown that the Angus model could provide an explanation of Hubble tension by enhancing cosmic variance, which also better accounts for the observed KBC void [53]. The discrepancy among the values of the Hubble constant obtained through the measurement of baryonic acoustic oscillations and the cosmic microwave background (67 km/s/Mpc) and the values deduced from the observations of cosmic candles (72–75 km/s/Mpc), also known as Hubble tension falsify $\Lambda$CDM at $7.09\sigma$, but considering these and several other observables leads to only $2.53\sigma$ for the MOND-hot dark matter hybrid cosmology of Angus [52,53]. We should also mention that this model also better accounts for the massive high-redshift interacting cluster “El Gordo”, which falsifies $\Lambda$CDM at $6.16\sigma$ confidence but could well have formed in MOND [54].

Among many other proposals using modified gravity to mimic the effects conventionally attributed to CDM, we can cite non-minimal gravitational couplings with matter [55] and bimetric theories [56].

A fully relativistic theory that has reaped some success in applications to astrophysics and cosmology is the scalar-tensor-vector gravity (STVG) or Modified Gravity (MOG) of J. Moffat and collaborators [57,58,59,60,61]. This theory is based on an action principle with an extra vector field and three scalar fields: the gravitational “constant”, the mass of the vector field and the coupling parameter of the Maxwell-Proca Lagrangian for the vector field. This theory accounts for the anomalous rotation curves of galaxies [58], the Tully-Fisher law [59] and the acoustic peaks of the Cosmic Microwave Background, accelerating expansion of the universe as well as the matter power spectrum [60]. All this is achieved without taking into account any form of dark matter [57]. Unfortunately, the fact that electromagnetic and gravity waves arrived simultaneously in the GW170817 event falsifies MOG theory in its present form [62]. The reason is that gravity waves propagate along ordinary geodesics in this theory but electromagnetic radiation follows a different modified trajectory as a result of the “emulated dark matter”. Recently, it has also been shown that the observed velocity dispersion profile of the ultra-diffuse galaxy Dragonfly 44 do not agree with the predictions of MOG at a $5.49\sigma$ level [63]. Another issue that contributes to the falsification of MOG is the fact that the dynamical discrepancies in galaxies set in below a certain acceleration—as in MOND—rather than beyond a fixed distance—as assumed in MOG with its fixed length scale [49].

From a phenomenological point of view, there have been some proposals to revise the validity of Newton’s gravity force on scales much larger than Solar System. As early as 1963, Finzi questioned the validity of Newton’s law at long distances. He proposed that it should be replaced by a force decreasing inversely with distance to the source instead of the inverse square [64]. Similar models were developed by Tohline [65], Kuhn & Kruglyak [66] and, more recently, by other authors [67,68,69].

From the previous fast review of the present status of the DM hypothesis and the modified gravity alternatives, we deduce that the situation in this research field is characterized by a lack of definitive convincing arguments in favour of any of the two concepts or particular theories. All the models have some successful predictions but also problematic comparisons with some observations and experiments.

To resolve this impasse it seems that new observations are necessary, especially on the scale of stars and planetary systems. It seems unlikely that a new theory of gravity could be established solely on the basis of predictions for galactic and cosmological scales. General Relativity gained acceptance through classical tests on the Solar System scale [70]—the same is expected to happen for an extended theory of gravity.

An observation of this kind may have been extracted for a selected set of wide binary stars using the Hipparcos and Gaia catalogues [71,72]. Hernández et al. have shown that a phenomenon similar to the flat rotation curves of galaxies is also present on the much smaller scale of binary stars separated by more than 7000 AU. Later on, McCulloch and Lucio gave an explanation of this perplexing phenomenon in terms of the theory of Quantized Inertia [73]. In this theory—as proposed by McCulloch in 2007—inertia generates as a consequence of the Unruh radiation of accelerating objects in the interplay of the Rindler horizon of the accelerating object and the Hubble horizon of the Universe. This model also yields good results for the anomalous rotation curves of galaxies [74]. However, these claimed successes rely on major mathematical flaws in calculation conducted within the framework of this theory, fixing which would lead to very substantial corrections to the numerical results [75].

We should notice that this distance scale coincides with:

$${\beta }_{\mathrm{MOND}}=\sqrt{{\displaystyle \frac{GM}{a_0}}},$$

(2)

where $GM$ is the mass constant of the Sun and $a_0$ is the critical acceleration of MOND. If confirmed, this phenomenon cannot be easily accommodated in the dark matter paradigm. Firstly, all forms of DM are uniformly distributed on these scales, so the DM in a 7 kAU sphere would exert a negligible amount of Newtonian gravity. This distribution is uniform because the velocity dispersion of CMD is of order 200 km/s, greatly exceeding Solar System velocities. Moreover, we can estimate the magnitude of the gravitational pull of DM in the Solar system in terms of the acceleration of a test particle at the edge of a uniform—spherically symmetric—mass distribution of density $\rho$ and radius r:

$${\displaystyle a_{\mathrm{DM}}=\frac{4\pi }{3}G\rho r.}$$

(3)

Taking into account the upper bound for the density of DM in the Solar system as obtained by de Salas et al. [76] in solar masses per cubic parsec, $\rho <0.01M_{\odot }/{pc}^3$, we obtain from Equation (3) that the effect of DM at 7000 AU is $a<2\times {10}^{-16}$ m/s${}^2$. This is six orders of magnitude below the Newtonian acceleration at that distance and, consequently, we can fairly say that DM would have a negligible impact on wide binaries.

The objective of this paper is to propose a modified Newtonian model that accounts for the flat rotation curves both for the scale of wide binary stars and galaxies. A modified scalar potential, satisfying a generalized Klein-Gordon equation, is considered to be the source of an extra force with a mass-dependent distance scale.

The paper is organized as follows: In Section 2 we discuss the extended Newtonian gravity model with an extra potential (MNG model on the following). This model is applied to the wide binary anomaly in Section 3. In this section we also derive the Tully-Fisher relation from the non-Newtonian model and we consider the radial velocity perturbations of Proxima Centauri as a possible test of MNG. Finally in Section 4 we calculate the contributions of the extra force to the advance of the planets’ perihelia, and discuss the possibility of detecting these small effects in new detailed planetary ephemerides. The paper ends with a discussion in Section 5 and some conclusions in Section 6.

2. Modified Newtonian Gravity

Newton’s theory can be understood as a scalar non-relativistic model of gravity based on a potential function. General Relativity also has a Newtonian limit for weak fields and low velocities of the bodies in comparison with the speed of light. Therefore, a modification of the Newtonian limit would also imply that General Relativity must be adapted and extended in some way.

Newton’s scalar potential (per unit mass) would be written as follows:

$${\displaystyle \varphi =-G\frac{M}{r},}$$

(4)

where G is Newton’s constant and M is the mass of the source within radius r. Following Finzi [64], Tohline [65] and other authors [66], we will have an extra force term at large distance from the source. This term would decay as $1/r$ with distance to the center of the—spherically symmetric—source mass. Consequently, we will have an extra potential $\chi \propto lnr$. With the guidance of this behaviour for large separations and the fact that anomalous accelerations within the Solar system are extremely small, according to the latest ephemeris [77], we have proposed an extra potential that satisfies a non-linear Klein-Gordon equation of the form:

$${\nabla }^2\chi =\alpha e^{-\gamma \chi },$$

(5)

where $\alpha$ and $\gamma$ are constants for a particular distribution of source mass. We should notice that Equation (5) is only valid for an isolated system since the parameter $\gamma$ vary between systems. Consequently, it cannot be applied to a system of interacting galaxies and we will use it only for particular source mass distributions. The study of many body systems would require well-defined field equations that we will not attempt to propose in this paper.

Extra potentials satisfying a Klein-Gordon equation with variable mass are also proposed in other models of modified gravity—such as MOG [57,58,59,60] in which a Maxwell-Proca Lagrangian with an interaction field is considered—and also in Quantum Field Theory [78]. In our case the exponential form of the self-interaction potential is chosen in such a way that it produces a $1/r$ dependence of the extra anomalous force at large distances as expected if we want to explain the flat rotation curves of galaxies. Our model has two parameters: $\alpha$ corresponding to the strength of the coupling of the anomalous potential and $\gamma$ that would be related with the distance scale of this potential.

In a spherically symmetric case—such as the field generated by a—approximately—spherical star—we can write Equation (5) as:

$${\displaystyle \frac{1}{r^2}\frac{d}{dr}\left[{\displaystyle r^2\frac{d\chi }{dr}}\right]=\alpha e^{-\gamma \chi }.}$$

(6)

This equation cannot be solved analytically, but we can study the limits $r\to 0$ and $r\to \infty$. For $r\to 0$ we propose a solution in the form of a harmonic potential:

$${\displaystyle \chi \simeq \frac{\mathcal{G}}{6}\frac{r^2}{{\beta }^3},}$$

(7)

where $\beta$ is a parameter with dimensions of distance and value much larger than the distance to the farthest planet, and $\mathcal{G}$ is a constant with the dimensions and order of magnitude of Newton’s gravitational constant times the mass of the stars. In the Solar System, $\chi \ll 1$ because $r\ll \beta$. By substituting the proposal in Equation (7) into Equation (6), we find that

$${\displaystyle \alpha =\frac{\mathcal{G}}{{\beta }^3}.}$$

(8)

Similarly, for $r\gg \beta$ we will have the potential $\chi \to ,\mathcal{G}/\beta ln(r/\beta )$. This is a solution of Equation (6) if the parameter $\gamma$ is chosen as follows:

$${\displaystyle \gamma =\frac{2\beta }{\mathcal{G}}.}$$

(9)

Equation (6) can also be integrated numerically for intermediate values of r. With this aim it is easier to work with non-dimensional quantities by performing the following redefinitions:

$$\begin{array}{ccc}\hfill r& \to & \rho \equiv r/\beta ,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \varphi & \to & \psi \equiv \beta \chi .\hfill \end{array}$$

(11)

Using these definitions, the equation for the scaled potential, $\psi$, becomes:

$${\displaystyle \frac{1}{{\rho }^2}\frac{d}{d\rho }\left[{\displaystyle {\rho }^2\frac{d\psi }{d\rho }}\right]=e^{-2\psi }.}$$

(12)

For a spherically symmetric and—approximately point-like—distribution of mass, we will have (from Equation (7)), the boundary conditions: $\psi \left(0\right)=0$, ${\psi }^{\prime }\left(0\right)=0$. The result for the numerical integration of Equation (12) is shown in Figure 1. Although analytic integral of this differential equation can be obtained, it is useful to consider the approximation:

$${\displaystyle \psi \left(\rho \right)={\int }_0^{\rho }du\frac{{tanh}^2\left(u\right)}{u},}$$

(13)

In Figure 1 we compare the numerical integration of Equation (12) with the approximation in Equation (13). The agreement is quantitatively and qualitatively good, so we will use in the following the extra potential in Equation (13). From this potential we derive the extra anomalous force (that adds to the classical Newtonian gravitational force):

$${\displaystyle \mathcal{F}\left(r\right)=-\frac{\mathcal{G}}{{\beta }^2}\frac{\partial \psi }{\partial \rho }=-\frac{\kappa GM}{\beta r}{tanh}^2\left({\displaystyle \frac{r}{\beta }}\right).}$$

(14)

Here we have made the substitution $\mathcal{G}\equiv \kappa GM$, with $\kappa$ being a numerical constant of order of unity, G is Newton’s constant and M is the mass of the point-like source.

To apply the model in Equation (14) we need a postulate for the length scale parameter, $\beta$. This was already suggested by Acedo [69] in connection with a similar modified gravity model for the rotational curves of galaxies. The expression for $\beta$ corresponding to a point source was postulated to take the form:

$$\beta =\eta \sqrt{M}.$$

(15)

As in Ref. [69] the parameter $\eta$ is a constant $\simeq 1\mathrm{m}/{\mathrm{kg}}^{1/2}$. Notice that a similar mass-dependent scale appears in STVG [57,58,59,60] and also in MOND as $r_{\mathrm{M}}=\sqrt{GM/a_0}$ [49].

3. Wide Binaries and the Tully-Fisher Relation

In 2018, Hernández, Jiménez and Allen selected a group of 280 wide binaries studied by the Hipparcos satellite that was, a highly accurate astrometric experiment in space [71]. These stars were chosen to avoid false alignments with a proportion of false positives $\simeq 10$ %. The objective of these researchers was to test the Newtonian prediction for relative velocities of binary systems that, for stars of equal mass, is given by:

$$\Delta V_N={\left({\displaystyle \frac{2GM}{r}}\right)}^{1/2},$$

(16)

where M is the (identical) mass of the stars and r their separation. Surprisingly, these authors found that for separations larger that $0.034$ parsecs (pc) or 7000 AU there is consistency with the predictions of Equation (16) within the errors bars. However, for $r>7000$ AU the relative velocities tend to reach a “plateau” of constant relative velocities irrespective of the separation among the components of the binary system. This finding was confirmed by another careful selection of binaries in the Gaia Data Release 2 (DR2) catalogue [72].

If established, by the results of future analyses, this discovery could be of the utmost importance for the physics of gravity. The reason is that it provides a replica for the anomalous rotation curves of galaxies and the anomalies in the applications of the virial theorem to galactic clusters, but on a much smaller scale. This would mean that the anomalies associated with the concept of dark matter are more profound and require a reformulation of General Relativity to take into account the breakdown of Kepler’s third law on large scales. A common explanation of the cluster, galactic and wide binary anomalies in the context of a DM model seems unlikely because DM tends to distribute uniformly on the scale of binary stars, so there is very little of it on the scale of a wide binary.

The obvious alternative approach to DM that could explain the wide binary observations is MOND [47,48]. The distance at which Hernández et al. [71,72] found that Kepler’s third law is not valid anymore corresponds to a stellar orbital acceleration $a_0\simeq 1.2\times {10}^{-10}$ m/s${}^2$. This is the same critical acceleration at which the effect of DM is noticed in galaxies [47]. Following McCulloch & Lucio [73] we will derive the prediction for the relative velocity of a pair of equal mass binary stars in the Solar neighbourhood of the Milky Way. The acceleration of the Sun around the Galaxy $a_g\simeq {10}^{-10}$ m/s${}^2$. Milgrom’s condition for a circular orbit of two equal mass stars is given as:

$${\displaystyle \frac{Ma}{1+a_0/(a+a_g)}=\frac{G{M}^2}{r^2},}$$

(17)

where $a+a_g$ is the total acceleration of the stars (the orbital acceleration of one star around the other plus the acceleration around the Galactic center). This is the so-called external field effect in MOND, which significantly affects the predictions. This effect has recently been detected by comparison of disk galaxy rotation curves in different external fields, which were previously estimated based on their environment [79]. Galaxies in more crowded environments tend to have a declining rotation curve at large distance from the galactic center, while the rotation curve remains flat in more isolated environments. Taking into account that the centrifugal acceleration would be $a=2v^2/r$, we obtain from Equation (17) a biquadratic equation in v that yields the following result for the relative velocity of the binary system ($\Delta v=2v$):

$$\Delta v_{\mathrm{MOND}}=\sqrt{{\displaystyle \frac{GM}{r}-a_{g}r+\sqrt{{\left({\displaystyle a_{g}r-\frac{GM}{r}}\right)}^2+4GM(a_g+a_0)}}}.$$

(18)

We will see that the external field effect makes the predictions of MOND not very different from that of Newtonian theory. We should notice that the predicted orbital velocity of wide binaries in the Solar neighbourhood separated by more than 7 kAU exceeds the Newtonian prediction by about 20% [80].

Our modified Newtonian model would correspond to a total acceleration of a component of the binary system given by:

$${\displaystyle a_{\mathrm{MNG}}=\frac{GM}{r^2}+\kappa \frac{GM}{\beta r}{tanh}^2\left({\displaystyle \frac{r}{\beta }}\right),}$$

(19)

and with a centrifugal acceleration of $2v^2/r$ we obtain the relative velocity of the binary system:

$$\Delta v_{\mathrm{MG}}=\sqrt{{\displaystyle \frac{2GM}{r}+\kappa \frac{2GM}{\beta }{tanh}^2\left({\displaystyle \frac{r}{\beta }}\right)}}.$$

(20)

We now assume a binary system with two stars of the same mass, this mass being the mass of the Sun, $M=1.989\times {10}^{30}$ kg. Then for $\eta =0.927$ m/kg${}^{1/2}$ as proposed in Ref. [69] we obtain $\beta \simeq 8737$ AU (see Equation (15)). We also estimate $\kappa =3.109$. With these values we obtain—from Equation (20)—the predictions for the relative velocity of the components of this system at various distances as shown in Figure 2. The agreement with the relative velocities obtained by Hernández et al. [72] is good, including a possible bump in the curve for distances in the range $0.01$–$0.1$ pc. The asymptotic velocity predicted by MNG in case of large separations is given by:

$$\Delta v_{\mathrm{MG}}=\sqrt{{\displaystyle \kappa \frac{2GM}{\beta }}}=\sqrt{{\displaystyle \frac{M}{M_{\odot }}}}0.796\mathrm{km}/sr\gg \beta ,$$

(21)

where $M_{\odot }$ is the mass of the Sun.

On the other hand, the predictions of Newtonian theory and MOND in Equations (16) and (18), respectively, are clearly below the observations’ error bars for $r>0.3$ pc. At this point we should emphasize that the wide binary data at distances greater than about 1 pc is not reliable for use in the wide binary test because such binaries would be rather susceptible to disruption by passing stars. This is quantified by the Jacobi radius, which is about 350 kAU—or $1.69$ pc—for two Solar mass stars in the Solar neighbourhood of our Galaxy [82]. As a result, wide binaries with a very large separation may be unbound. For this reason, Banik and Zhao [80] recommended to only use data out to a separation of 20 kAU—$\simeq 0.1$ pc—and, certainly, data for separations larger than 1 pc would be contaminated by many false binaries. A more careful analysis including outlier rejection has been carried out by Banik [83]. In this study, the author concludes that the widest sky-projected separation bins should not be used. The problem is that at lower separations the observations do not distinguish MOND from Newtonian dynamics (or the MNG model in this paper). Nonetheless, the results of Pittordis and Sutherland [81] strongly rejects models where local wide binaries have a constant asymptotic orbital velocity and leave only MOND with external field effect as a valid alternative to DM on these scales.

3.1. The Tully-Fisher Relation

The Tully-Fisher (TF) relation is a phenomenological law—obtained from statistical analyses of the rotation curves and visible mass content of many galaxies. In 1977, the astronomers R. B. Tully and J. R. Fisher [50] published the result of their analysis showing that there is a relation among the maximum rotation velocity—achieved in the rotation curve of a spiral galaxy, and its mass content or luminosity [50,84]. This relation takes the following form:

$$M=q{V}_{\max }^{\nu },$$

(22)

where q and $\nu$ are constants. We assume that M is measured in solar masses and $V_{\max }$ in km/s. As a measure of the mass of the galaxy, Tully and Fisher used the optical luminosity which is directly related to the total stellar mass. This yields the so-called stellar TF relation [85]. Later on, some authors considered also the gaseous content of the galaxy (mainly in the form of hydrogen clouds detected through the 21 cm line of neutral hydrogen and H$\alpha$ spectra). If both the stellar mass and the gas content are taken into account in the evaluation of the total mass M, the baryonic TF relation is obtained. Torres-Flores et al. [84] have given the following values for the parameters in this baryonic TF relation:

$$\begin{array}{ccc}\hfill q& =& {10}^{2.21\pm 0.61},\hfill \\ \hfill \nu & =& 3.64\pm 0.28.\hfill \end{array}$$

(23)

The TF relation encapsulates the anomaly found in the rotation curves of spiral galaxies by predicting the asymptotic constant rotation velocity of stars and gas far off the galactic centers in terms of the total mass, or vice versa. The corresponding empirical relationship for elliptical galaxies is known as the Faber-Jackson relation [86] , so that these anomalies are found in all kind of galaxies.

A well-known consequence of MOND is that it predicts a TF relation with exponent $\nu =4$ [47]. This is largely considered to be an early success of the MOND approach. Contrarily to what is usually thought, modified Newtonian gravity can also predict a TF relation if it incorporates a mass-dependent length scale [69]. For example, for a mass of the Milky-Way estimated as $75\times {10}^9$ solar masses we obtain $\beta =11.6$ kpc. At distances much larger than this critical value, and assuming that mass is concentrated at the core of the galaxy, we can approximate the anomalous acceleration in our model (deduced from Equation (19)) as:

$${\displaystyle a(r\gg \beta )\simeq \kappa \frac{GM}{\beta r},}$$

(24)

where M is the total mass of the galaxy (stars and gas clouds) and $\beta$ is proportional to the square root of this mass as given in Equation (15). By equating this asymptotic acceleration to the centrifugal one we obtain the relation:

$$M={\left({\displaystyle \frac{\kappa G}{\eta }}\right)}^2{V}^4.$$

(25)

This corresponds to TF relation with exponent $\nu =4$ and coefficient $q={(\kappa G/\eta )}^2$. For $\kappa =3.109$ and $\eta =0.927$ m/kg${}^{1/2}$ this coefficient would take the value $q\simeq 5.01\times {10}^{-20}$ (if M in Equation (25) is measured in solar masses and V in km/s). From Equation (23) the fitted value of the proportionality constant in the TF relation would be in the range $7.61\times {10}^{-22}<q<1.26\times {10}^{-20}$, although there is large uncertainty. We conclude that the modified Newtonian gravity model discussed in this paper can explain both the anomalous relative orbital velocity of wide binaries and the TF relation. It would therefore would bridge the gap among the gravitational anomalies ranging from stellar to galactic scales.

3.2. Radial Velocity Perturbations of Proxima Centauri

Proxima Centauri is a red dwarf star that is a member of the Alpha Centauri system. This red dwarf is also the closest star to the Sun. It is currently at a distance of 12950 AU from the Alpha Centauri A and B components. For this reason, it is also an ideal candidate for the study of the wide binary anomaly. In 2018, Banik and Zhao [80] proposed an analysis of the orbit of Proxima Centauri throughout a period of 5 years using the future Theia Space Telescope [87]. It was shown that if this new mission achieves a precision of 1 $\mathsf{\mu }$ as during that period of time it could resolve the sky-projected velocity components of these star with enough accuracy to elucidate whether there is an anomalous trajectory—of the kind predicted by MOND, MNG or other models—or not. Astronomers have already been able to fit the Keplerian osculating orbit of Proxima’s Centauri at present [88]. The orbital parameters are listed on

Table 1. Orbital parameters for the osculating keplerian orbit of Proxima Centauri at present as given by Kervella et al. [88]. The angles are written in sexagesimal degrees.

Parameter	Value and Units	Plus Error	Minus Error
a Semi-major axis	8700 AU	$+700$ AU	$-400$ AU
e Orbital eccentricity	$0.5$	$+0.08$	$-0.09$
P Period	$547,000$ years	$+66,000$ years	$-40,000$ years
$T_0$ Epoch of periastron	$283,000$ years	$59,000$ years	$-41,000$ years
I Orbital inclination	${107.6}^{\circ }$	$+{1.8}^{\circ }$	$-{2.0}^{\circ }$
$\Omega$ Longitude asc. node	${126}^{\circ }$	$+5^{\circ }$	$-5^{\circ }$
$\omega$ Argument of periastron	${72.3}^{\circ }$	${8.7}^{\circ }$	${6.6}^{\circ }$
$RV$ Radial velocity	$-21.7$ km/s	$+0.027$ km/s	$-0.027$ km/s

.

In 2017, Iorio also suggested to use the radial velocity (RV) of a star as the key observable obtained by Doppler spectroscopy. This way the effects of General Relativity could also be detected in S2 (a star revolving around the massive black-hole at the center of the Galaxy) [89]. In this section we will apply Iorio’s method to the problem of detecting deviations from a Keplerian orbit in a wide binary as a consequence of nonstandard perturbations as the one predicted by our MNG model. The ideal Keplerian orbit of Proxima Centauri would be given by:

$${\displaystyle r\left(f\right)=\frac{a(1-e^2)}{1+ecosf},}$$

(26)

where f is the true anomaly and $r\left(f\right)$ is the distance to the center of mass of Alpha Centauri. The relation among time and the true anomaly is also useful:

$${\displaystyle \frac{dt}{df}=\frac{{(1-e^2)}^{3/2}}{n_b{(1+ecosf)}^2},}$$

(27)

where $n_b=\sqrt{\mu a^{-3}}$ is the keplerian mean motion (with $\mu$ being the total mass parameter of the binary). The mass of Proxima Centauri is estimated as $M_{\mathrm{Proxima}}=0.1221M_{\odot }$ and the total mass of Alpha Centauri A and B is $M_{\mathrm{total}}=M_A+M_B=2.007M_{\odot }$.

The perturbing force corresponding to the MNG model in this paper is according to Equation (14):

$${\displaystyle A_{\rho }=-\kappa \frac{\mu }{\beta r}{tanh}^2\left({\displaystyle \frac{r}{\beta }}\right),}$$

(28)

with $\kappa =3.109$ and $\beta =12,377.6$ AU, as deduced from Equation (15). Following Iorio [89] we find that the nonzero contributions to the perturbation of the orbital elements are the ones corresponding to the semi-major axis, the orbital eccentricity and the argument of the periapsis:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{da}{dt}}& =& {\displaystyle \frac{2e}{n_b\sqrt{1-e^2}}{A}_{\rho }sinf,}\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{de}{dt}}& =& {\displaystyle \frac{\sqrt{1-e^2}}{n_{b}a}{A}_{\rho }sinf,}\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d\omega }{dt}}& =& {\displaystyle -\frac{\sqrt{1-e^2}}{n_{b}ae}{A}_{\rho }cosf.}\hfill \end{array}$$

(31)

For $\kappa =a$, e, $\omega$ the total variation of the element in the interval corresponding to a variation of the true anomaly from $f_0$ to f is given by:

$${\displaystyle \Delta \kappa \left(f\right)={\int }_{f_0}^{f}d{f}^{\prime }\frac{dt}{d{f}^{\prime }}\frac{d\kappa }{dt},}$$

(32)

where $dt/d{f}^{\prime }$ is evaluated in Equation (27). We also need the perturbation of the mean anomaly as follows:

$${\displaystyle \Delta \mathcal{M}\left(f\right)=\Delta ϵ\left(f\right)-\Delta \overline{\omega }-\frac{3n_b}{2a}{\int }_{f_0}^{f}d{f}^{\prime }\Delta a\left(f^{\prime }\right)\frac{dt}{d{f}^{\prime }},}$$

(33)

where $\Delta ϵ$ is the perturbation of the mean anomaly at epoch and $\Delta \overline{\omega }$ is the perturbation of the longitude of the pericenter calculated as follows:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d\overline{\omega }}{dt}}& =& {\displaystyle -\frac{\sqrt{1-e^2}}{n_{b}ae}{A}_{\rho }cosf,}\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dϵ}{dt}}& =& {\displaystyle \frac{e^2}{1+\sqrt{1-e^2}}\frac{d\overline{\omega }}{dt}-\frac{2}{n_{b}a}{A}_{\rho }\left({\displaystyle \frac{r}{a}}\right),}\hfill \end{array}$$

(35)

the finite perturbations being calculated also with Equation (32). Finally the perturbation of the RV for Proxima Centauri would be given by:

$$\begin{array}{ccc}\hfill \Delta RV\left(f\right)& =& {\displaystyle \frac{M_{\mathrm{Proxima}}}{M_{\mathrm{total}}}\left\{{\displaystyle -\frac{n_{b}sinI(ecos\omega +cos(\omega +f)}{2\sqrt{1-e^2}}\Delta a\left(f\right)}\right.}\hfill \\ \hfill & +& n_{b}asinI\left(4cos(2f+\omega )-ecos(f-\omega )+4ecos(f+\omega )\right.\hfill \\ \hfill & +& \left.ecos(3f+\omega )\right)/4{(1-e^2)}^{3/2}\Delta e\left(f\right)\hfill \\ \hfill & -& {\displaystyle \frac{n_{b}asinI(esin\omega +sin(\omega +f)}{\sqrt{1-e^2}}\Delta \omega \left(f\right)}\hfill \\ \hfill & -& \left.{\displaystyle \frac{n_{b}a{(1+ecosf)}^{2}sinIsin(\omega +f)}{{(1-e^2)}^2}\Delta \mathcal{M}\left(f\right)}\right\}.\hfill \end{array}$$

(36)

Taking into account the epoch of periastron in Table 1 we can estimate the true anomaly today as $f_0=3.0995665$ by using the standard equations connecting the true anomaly with the eccentric anomaly and time [90]. After ten years the true anomaly becomes $f=3.0996107$. By using the Equations (26)–(36) with the parameters in Table 1 and considering a interval of time of ten years from now we obtain the results in Figure 3 for the perturbation of the radial velocity as a consequence of the non-Newtonian force in the MNG model.

We notice that the magnitude of this perturbations is below the error bars for the determination of Proxima Centauri’s radial velocities (as shown in Table 1). Consequently, it would not be possible to test the MNG model with Proxima Centauri in the near future using radial velocities alone.

4. Modified Newtonian Gravity and the Anomalous Perihelion Precessions in the Solar System

Modern ephemerides take into account a vast amount of data of increasing accuracy. Radar, laser ranging and the tracking of interplanetary spacecraft is constantly improving these ephemerides. Pitjeva et al. [91] already developed in 2008 an ephemerides model (EPM2008) that took into account a dataset of planetary positions for the period 1913–2007 including also the perturbation by trans-Neptunian objects, the solar quadrupole, the 301 biggest asteroids and a ring that simulates the gravitational effect of the mass of the smaller asteroids. This ephemeris also incorporates the corrections predicted by the Post-Newtonian approximation of General Relativity.

This model, and others developed later, computed the perihelion precession and any deviation from the predictions of General Relativity [77,91,92,93,94]. Although any anomalies beyond General Relativity would be very small (as the observations are telling us), there is a chance of discovering of a minute extra perihelion advance (or regression) in the range of a few milliarcseconds per century. In particular, for EPM2008 an anomalous retrograde motion of the perihelion of Saturn was found. This anomaly was estimated as $-6\pm 2$ milliarcseconds per century, but it was not confirmed by other ephemerides [94]. Notwithstanding the still large uncertainties, we cannot dismiss that a discrepancy with General Relativity could be discovered in the near future.

For a radial perturbing force, $\mathcal{R}$, the argument of the pericenter (as measured in the orbital plane from the line of nodes) would change as follows [95]:

$${\displaystyle \frac{d\omega }{dt}=-\sqrt{{\displaystyle \frac{p}{GM}}}\frac{1}{ϵ}\mathcal{R}cosf,}$$

(37)

where $p=a(1-{ϵ}^2)$ is the semi-latus rectum (a being the semi-major axis of the orbit and $ϵ$ the orbital eccentricity) and f is the true anomaly measured from pericenter. For an ideal keplearian orbit there is a relation among time and the true anomaly [96] as follows:

$${\displaystyle dt=\frac{T}{2\pi }\frac{{(1-{ϵ}^2)}^{3/2}}{{(1+ϵcosf)}^2}df,}$$

(38)

where T is the orbital period given by Kepler’s third law:

$$T=2\pi a^{3/2}{\left(GM\right)}^{-1/2}.$$

(39)

From Equation (14) we now have that the non-Newtonian extra force can be written as:

$${\displaystyle \mathcal{R}=-\frac{\kappa GM}{{\beta }^3}r,}$$

(40)

where M is the mass of the Sun and we have taken into account that for the planets $r\ll \beta$. We should also consider that for Keplerian orbits, $r=p/(1+ϵcosf)$. This would be useful to estimate the extra force [95].

From Equations (37), (38) and (40) we can obtain the perihelion advance in terms of the true anomaly:

$${\displaystyle \frac{d\omega }{df}=\frac{T}{2\pi ϵ}\frac{\kappa }{{\beta }^3}{p}^{3/2}{\left(GM\right)}^{1/2}cosf\left(1-3ϵcosf+\mathcal{O}\left({ϵ}^2\right)\right).}$$

(41)

As is usual in perturbation theory, we now calculate the average over a complete revolution $0\le f\le 2\pi$ taking into account that $\langle cosf\rangle =0$ and $\langle {cos}^{2}f\rangle =1/2$. We obtain the following result for the anomalous perihelion retrograde motion (per revolution):

$${\displaystyle \Delta \omega =-\frac{3}{2}\kappa {\left({\displaystyle \frac{a}{\beta }}\right)}^3,}$$

(42)

where we have used Kepler’s third law in Equation (39) to simplify the result. For the Solar System we proposed in the previous Section that $\kappa =3.109$ and $\beta =8737$ AU. From the values of the orbital period and semi-major axes of the planetary orbits [97], we get an anomalous advance of the perihelion by $-0.034$ mas per century for Mercury, $-0.271$ mas cy${}^{-1}$ for Mars, $-1.713$ mas cy${}^{-1}$ for Jupiter and $-4.307$ mas cy${}^{-1}$ for Saturn. This is to be compared with some ephemerides such as EPM2008 [91], EPM2011 [93] or the most recent EPM2017 as analyzed by Iorio [98]. The comparison is shown in

Table 2. Prediction of the MNG model and MOND for the extra anomalous advance of the perihelion of the planets in comparison with the observations in the EPM2008 and EPM2011 ephemerides. In the case of EPM2017 only the formal uncertainties are given as deduced by Iorio  [98].

Planet	MNG (mas cy${}^{-1}$)	MOND	EPM2008	EPM2011	EPM2017
Mercury	$-0.034$	$0.02$	$-3.6\pm 5.0$	$-2\pm 3$	$0.008$
Mars	$-0.271$	$-0.09$	$0.1\pm 0.5$	$-0.02\pm 0.037$	$0.003$
Jupiter	$-1.713$	$-0.65$	−	$58.7\pm 28.3$	$33.9$
Saturn	$-4.307$	$3.12$	$-6\pm 2$	$-0.32\pm 0.47$	$0.067$

. The predictions of MOND—using the standard MOND function $\mu \left(y\right)=y/\sqrt{1+y^2}$—is also given [99].

There has been some interest about the significant measurement of an anomalous perihelion precession for Saturn in EPM2008 [94], although this was not confirmed by other ephemerides. On the other hand, We must notice that EMP2011 disagrees from our predictions for Mars for more than $6\sigma$ and from that of Saturn for more than $8\sigma$. Disagreement is even larger with the smaller uncertainties of EPM2017. A similar discrepancy is found when the anomalous quadrupolar corrections predicted by MOND—along the direction of the external galactic field—are taken into account [99].

Therefore, even if our model predicts such an anomaly as found in EPM2008, both in sign and order of magnitude, it is premature to obtain any conclusion from this result because the uncertainties in the ephemerides are large and change from one model to another. However, if the accuracy of these models increase, we may have a confirmation of an anomaly in the precession of the planetary orbits that, as shown in this paper, could be connected with the breakdown of Kepler’s third law in wide binary systems and in the outskirts of galaxies. The accuracy of ephemerides is increasing continuously thanks to the data recollection of high accuracy measurements using radio and laser techniques. In particular, Lunar laser ranging—which uses the timing of reflected laser beams on retroreflectors placed on the surface of the Moon by the Apollo missions [100]—is providing data of unprecedented accuracy in celestial mechanics.

5. Discussion

The gravitational anomaly discovered by Fritz Zwicky almost a century ago has presently an unquestioned importance in our current understanding of the Universe [1,5,15]. A different issue concerns the explanation and role of this anomaly in astrophysics and cosmology. Most researchers adhere to Zwicky’s initial hunch that the failure of the virial theorem in clusters was caused by a large mass of dark matter that was not taken into account in the calculations. The scientific community kept faithful to Zwicky’s idea when a similar anomaly was discovered in galaxies by Rubin et al. [3] and Faber and Gallagher [86]. Since then the concept of DM has been incorporated in the $\Lambda$CDM model [5,6,7,8,9] where it provides some explanatory power in connection with the age of the Universe, the spectrum of fluctuations of the Cosmic Microwave Background and the formation of structure.

Despite this success, the DM idea has not convinced all researchers in the field and alternative explanations have been keenly investigated for many years. As early as 1963, Finzi suggested that Newton’s law of gravity could not be valid for galactic scales. In 1983, Milgrom proposed his—much debated—MOND model to explain the onset of the anomaly in galaxies below a critical orbital acceleration [47,48,49]. Later on, Moffat and collaborators produced a whole relativistic theory of modified gravity that could address many problems in cosmology [57,58,59,60].

The interest in these modified gravity theories is not merely academic but is encouraged by the fact that dark matter cannot always explain the observations. The most basic problem of the DM idea is the lack of experimental evidence of any of the many possible candidate particles that could be the main component of DM [15]. The DAMA/LIBRA and similar experiments are large underground particle detectors formed by sodium iodide scintillation detectors that look for weakly interacting massive particles in the Galactic halo [16,17]. The key signature of WIMPs is expected to be an annual modulation of the detection signal that has indeed been found, but could have a conventional explanation [18]. On the other hand, other experiments have not been able to reproduce the DAMA/LIBRA and DAMA/NaI results [19,20]. Similarly, for the sterile neutrino proposal the possibility of a decay channel into X-rays has been suggested [41,42] and a recently discovered $3.5$ keV feature in the X-ray spectra of some galaxies has been associated with this decay. A lot of wishful thinking is connected with this premature identification [43,44,45]. Other issues with DM that are not satisfactorily solved within the context of the $\Lambda$CDM model, are the cuspy-halo problem [21] and the missing satellites problem [23,24,25]. Equally important is the satellites plane problem [26] consisting on a very unlikely distribution of satellite galaxies in the Milky Way, Andromeda galaxy and other if the $\Lambda$CDM were correct.

Although neither the DM nor the modified gravity adherents can claim a complete solution to the whole riddle posed by the gravitational anomalies, it would be of great interest if any observation could elucidate the solution by giving strong evidence in favour of one of the two concepts. It seems that such a crucial experiment could have already been done in the form of the highly accurate astrometric catalogues provided by the Hipparcos and Gaia missions [71,72]. By using a careful selection of wide binary stars in these catalogues, Hernández et al. [71,72] have found that the relative orbital velocity of binary stars separated by $>7000$ AU violates Kepler’s third law within the margin of the error bars. This observation fits well in the history of the gravitational anomalies that gave rise to the concept of dark matter: In the first place, Zwicky discovered that galaxies in galactic clusters move too fast for the total luminous matter content [1]. Fifty years later it was found that a similar anomaly is present in galaxies [3,4], and after another forty years the same anomaly emerges in the context of binary stars [71,72]. If this is the case, it seems more natural to attribute the anomalies to a modification of gravity on large scales instead of DM because this breakdown of Kepler’s law appears to be universal and is found on very different scales. Both cold dark matter (WIMPS) and warm dark matter in the form of sterile neutrinos should not aggregate around stars to explain the wide binaries anomaly. Other forms of DM that have been discussed in the literature, could perhaps help in understanding this new anomaly such as dissipative dark matter [101]. However, attributing new properties to DM every time a new observation does not fit into the model is a poor scientific procedure.

6. Conclusions

In this paper we have studied a model of gravity with two scalar fields: one of them being the classical Newtonian potential and the other one appearing as the solution of a generalized Klein-Gordon equation. The second potential gives rise to a force that decays as—$\kappa G/\left(\beta r\right)$ with the separation r of the particles (for $r\gg \beta$). The scale of this force, $\beta$, depends on the mass of the source of the field and $\kappa G$ is the gravitational constant for this new force. The main thesis of this paper is that by taking $\beta =\eta \sqrt{M}$, where $\eta$ is a new constant of order $1\mathrm{m}{\mathrm{kg}}^{-1/2}$ and M is the mass of the source, we can explain both the wide binaries anomaly and the Tully-Fisher relation [50,84] that encapsulates the general features of the rotation curves of galaxies. Another key question concerning proposals of a modification of General Relativity is the fact that the tests of this theory in the Solar System are very precise [70]. Any proposal of this kind should agree with the bounds on perturbations provided by the latest ephemerides [91,92,93]. Of course, it would be more interesting if these ephemerides disclose a new anomaly that could be tested within the context of the new theory. The model discussed in this paper predicts an extra anomalous precession of the planetary orbits that could be proved right in the future if the accuracy of these ephemerides is increased to a fraction of milliarcsecond per century. However, we should emphasize that the latest ephemerides EPM2017 [98] do not confirm earliest findings and the MNG model discussed in this paper is strongly disfavoured if the new results are established.

Thanks to the observations of increasing accuracy that are being obtained in astrophysics and cosmology, we are now in a better position to propose a fully relativistic theory with General Relativity as its natural limit in the regime of low velocities, weak fields and small distances scales in comparison with the size of the Solar system. This theory should agree with the propagation of electromagnetic and gravitational waves along the same trajectories [62]. Unlike the standard theory it is expected to explain both the wide binary anomaly [73] and the anomalous rotation curves of galaxies [49]. We must also observe that the wide binary anomaly is far from being a solid fact in astronomy. Recently, Pittordis and Sutherland [81] have performed simulations for a large number of possible binary systems in MOND and Newtonian gravity and their results strongly reject any model—such as MOND without external field effect or MNG—that predicts an asymptotic constant orbital velocity for wide binary stars at large separations. In this context, MOND with external field effect would still be the best model of modified gravity. On the other hand, we must take into consideration its problems in the case of galactic clusters [52] as well as the large extra planetary precessions predicted for the Solar system [99].

Despite its shortcomings, the model discussed in this paper could provide a starting point to a theory of this kind and give support to the hypothesis of modified gravity over the traditional concept of dark matter. A debate that should be sorted out by means of the interplay between theory and observations instead of clinging to a concept that may have been useful but wrong.

An important next step would be to calculate the predicted velocity distribution of wide binaries in different gravity theories as a function of separation and mass, which allows for a much more detailed comparison with observations than considering only the velocity dispersion. The work of Pittordis & Sutherland [81] has the obvious implication that Monte Carlo—or similar methods—should be used to determine the expected distribution of the observable parameter $v_p/v_N\left(r_p\right)$, where $v_p$ is the projected velocity and $v_N\left(r_p\right)$ is the Newtonian prediction, i. e., $v_N\left(r_p\right)=\sqrt{GM/r_p}$ for wide binaries of total mass M and sky-projected separation $r_p$. This will clarify whether the peak of the distribution remains close to 1 even for a sky-projected separation $r_p=20$ kAU, as would be needed to match observations. Work along this line is currently under progress and it will be published elsewhere.