Absence of a fundamental acceleration scale in galaxies
Abstract
The Radial Acceleration Relation[1] confirms that a nontrivial acceleration scale can be found in the average internal dynamics of galaxies. The existence of such a scale is not obvious as far as the standard cosmological model is concerned, and it has been interpreted as a possible sign of modified gravity[2, 3]. The implications could be profound: it could in principle explain galactic dynamics without large amounts of yetundetected dark matter[4, 5] and address issues that the standard cosmological model faces at galactic scales[6, 7, 8, 9, 10]. Here, we consider 193 disk galaxies from the SPARC[11] and THINGS[12, 13] databases and, using Bayesian inference, we show that the probability of existence of a fundamental acceleration that is common to all the galaxies is essentially zero: the value is smaller than or, equivalently, the null hypothesis is rejected at more than 10. We conclude that the acceleration scale unveiled by the Radial Acceleration Relation is of emergent nature, possibly caused by a complex interplay between baryons and dark matter. In particular, the MOND theory[14, 15, 16, 17], or any other theory that behaves like it at galactic scales, is ruled out as a fundamental theory for galaxies at more than 10.

Center for Astrophysics and Cosmology, CCE, Federal University of Espírito Santo, 29075910, Vitória, ES, Brazil.

Department of Physics, CCE, Federal University of Espírito Santo, 29075910, Vitória, ES, Brazil.

Dipartamento di Fisica e Astronomia, Università di Catania, Viale Andrea Doria 6, 95125 Catania, Italy.

Institute of Modern Physics, Chinese Academy of Sciences, Post Office Box 31, Lanzhou 730000, People’s Republic of China.

INFN sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy.

Department of Physics, Bu Ali Sina University, Hamedan, Iran.
This version does not match the published version at Nature Astronomy.
The published version can be freely found at https://rdcu.be/ZXNT.
Dark matter is currently one of the main mysteries of the universe. There are many strong indirect evidences that support its existence, from galactic rotation curves and galaxy cluster dynamics to cosmological structure formation and cosmic microwave background anisotropies, but there is yet no sign of a direct detection[4, 5]. Moreover, at the scales of galaxies, there is tension between the theoretically expected dark matter distribution in the universe (from the standard cosmological model, CDM) and its indirectly observed distribution[6, 7, 8, 9, 10]. Therefore, phenomena associated to dark matter have a chance of serving as a window towards new physics.
Among the observational correlations whose explanation may be simpler within alternative theories of gravity, the Radial Acceleration Relation[1] (RAR) stands out[2]. It is a sharp correlation between the observed acceleration along the galactic radius (found from the redshift of atomic or ionised hydrogen) and the Newtonian acceleration of baryonic matter (found from the mass profiles that are derived from the observed surface brightness of the different galaxy components). Both these accelerations can be derived from observational data independently of any hypothesis on dark matter. From the CDM perspective, recent results using complex largescale hydrodynamical simulations[18, 19] indicate that CDM may reproduce such correlation. This work supports this line of investigation, but we stress that, within the standard model, the RAR is far from obvious. The RAR is also known as “Mass DiscrepancyAcceleration Relation” (MDAR).
The Modified Newtonian Dynamics (MOND)[14, 15, 16, 17] – although cannot explain so much of the observed universe as the CDM model does – correctly predicts the existence of certain features that are not so clear within the standard model, such as the TullyFisher relation[20, 21, 1]. An essential assumption of this approach is the existence of a fundamental acceleration scale, given by a constant , such that the gravitational force decays with distance slower than Newtonian for accelerations about or smaller than . More precisely, MOND was developed based on the following correction to the Newtonian acceleration in the context of disk galaxies[14, 15, 16]:
(1) 
where is the Newtonian acceleration, is the one expected to be the physical acceleration and is called the interpolating function. This function is such that for large accelerations () one would find Newtonian gravity (), while for small ones one gets . There is not much room to change these limiting cases since is required by Solar System dynamics and is essentially demanded by the TullyFisher relation in the absence of dark matter[22, 23]. The function controls the smoothness of the transition between Newtonian gravity and the regime of very small accelerations. A too sharp transition between these regimes will lead to issues on galaxy dynamics, since observational data indicates that the transition is smooth along the galaxy radius. While a too smooth transition is more prone to be observed via local tests of gravity[16]. The Simple Interpolating Function faces difficulties with Solar System constraints[24], while the Standard Interpolating Function has no issues with the Solar System, but performs worse with galaxies[25, 26]. These two interpolating functions are the most used ones. The RAR suggests another interpolating function, which comes directly from its results and we call here RARinspired Interpolating Function. The latter is an intermediate case for the range of Newtonian accelerations probed by the rotation curves of galaxies. In this work, we consider these three interpolating functions (see Methods for their expressions). As we will discuss, our conclusions are robust as far as these three interpolating functions are concerned, and likely as far as any other interpolating function that is compatible with the RAR.
The standard procedure to measure from galaxy rotation curves is as follows[26, 16]: first one performs fits on several galaxies considering as a free parameter, then one estimates the true value by taking the median (an estimator more robust than the mean). The quoted error is the error on the median[26], which scales as where is the size of the galaxy sample. The result of this procedure applied to the SPARC sample is presented in the Supplementary Material. Proceeding in this way, one will always find a best value for each interpolating function, but the existence of a fundamental universal is not tested, it is assumed from the beginning.
MOND is commonly criticized at larger than galactic scales, where its predictions fail to hold[27], or some form of additional matter seems to be necessary[28, 29]. There is hope that a more complete theory, that becomes similar to MOND at galactic scales, may have better results at all scales[30, 31, 32, 3]. For the Milky Way alone, assuming the Standard Interpolating Function and a large class of viable baryonic models, it was shown that the inferred in the Milky Way is incompatible with the one inferred from external galaxies at more than[33] . This alone is an interesting result, but it has relevant limitations since this conclusion does not hold for the Simple Interpolating Function and it depends on the analysis of a single galaxy (a special galaxy, but one nonetheless). There are other works that criticize MOND by presenting data that suggest a non constant value, but these criticisms were not based on the Bayesian analysis of a large sample of galaxies[34] (see also references therein).
Here, we test whether galaxy data is consistent with a fundamental acceleration scale. If it is not, then the RAR is an emergent phenomenon that appears when diverse galaxies are stacked together. In this case, MOND, as any other theory that behaves like it at galactic scales, is inconsistent with the data.
A straightforward way to assess the existence of a fundamental is to perform a full Bayesian analysis and obtain the posterior probability distribution of for each galaxy in the sample. One can then test if the various posteriors are compatible with each other. Bayesian inference has the great advantage that parameters are never fixed to some central or bestfit values and the posterior on is obtained by marginalizing over the other parameters. As such, we do not bias the result in any way: any possible degeneracy, expected or unexpected, linear or not, will be correctly taken into account. Furthermore, we do not linearize the likelihood with respect to the parameters and we obtain the confidence intervals by integrating the posterior distribution without adopting any Gaussian approximations. Further details of the analysis (including a recap of Bayesian inference) can be found in the Supplementary Material.
In order to obtain the (unmarginalized) posterior on , stellar masstolight ratios and (for the disk and the bulge, respectively) and galaxy distance , we adopt the following flat priors: , constrained to lie within of the observed distance, and and constrained to lie within a factor of their expected values. Changes on distance and masstolight ratios in general have significant impact on MOND results, since its dark matterlike effects come directly from the usual matter distribution. Hence, to properly evaluate MOND, it is important to consider a wideenough prior on these parameters. Our choice is of the same order, or larger, than the expected ranges for the SPARC data[1]. It is important to stress that stellar masstolight ratios and galaxy distance are treated here as nuisance parameters that model the impact of possible systematics. More precisely, as shown in the Supplementary Material, stellar masstolight ratios and galaxy distance are correlated with each other and with . Therefore, by marginalizing over the these nuisance parameters we are effectively increasing the error bars in a manner that is proportional to the widths of the priors we adopt.
Our main results are shown in Figure 1 and Table 1. Figure 1 shows the posterior distributions on , including their maxima and the 1, 3, and 5 confidence ranges (credible intervals). These results use no Gaussian or linear approximations. All the 175 SPARC galaxies[11] were analysed, but only 100 of them passed the quality cuts (see Methods) and were considered for the analyses of Figure 1 and Table 1. As it is evident, there is no value of that is inside all the intervals. Figure 1 also shows the global best fit of (see Methods). It implies that 13 galaxies, among the selected 100 galaxies, are incompatible with the global best fit of by more than . This issue is further detailed in Table 1. The existence of a fundamental universal acceleration scale can be discarded with a confidence of 10 or more, as detailed in Methods.
The results above are quite robust. As described in Methods, our conclusions remain unchanged after applying strict quality cuts (we removed 55% of the sample with the most strict quality criteria, see the Supplementary Material for the plot). Only if there were a strong and systematic underestimation of the errors then the results above could be somehow questioned. However, such a large increase would raise several other issues on rotation curve analyses. Among others, the RAR analysis itself indicates that the average expected errors are about or larger than the observed dispersion[1], hence increasing the errors would not improve the RAR significance. Moreover, in order to address possible systematic issues that correlate with galaxy luminosity or surface brightness, we also evaluate correlations of our Bayesian results with these parameters, as detailed in the Methods section and in the Supplementary Material. We find no significative correlation capable of changing our results.
The results shown in Figure 1 were found using the RARinspired Interpolating Function, but the Simple and the Standard interpolating functions yield similar results, as summarized in Table 1 and explicitly shown in the Supplementary Material. This is due to the fact that, as shown in Figure 2, the three interpolating functions are inside the dispersion of the RAR. Other possible interpolating functions cannot differ too much from the three interpolating functions we consider, otherwise they will not satisfy the RAR. In the Supplementary Material we also perform the same analysis for some galaxies of the THINGS sample[13]. Albeit with less galaxies, the analysis leads to the same conclusion: there is no common acceleration for all the galaxies.
Interpolating  Bayesian approach  

function  incompatible galaxies (expected)  
()  ()  ()  
Standard  12.902  67%  30%  13% 
Simple  12.918  66%  28%  11% 
RARinspired  12.959  64%  28%  13% 
Table 1 summarizes the visual results of Figure 1 by presenting the percentage of galaxies that are incompatible with the global best fit for at different levels. Reference values are given in parenthesis: there are clearly many more incompatible galaxies than expected.
Concluding, the RAR indicates the existence of an acceleration scale that emerges from the analysis of several galaxies, when their data is stacked together. On the other hand, from the combination of the individual galaxy analyses, we show, using Bayesian inference, that this emergent acceleration cannot be considered a fundamental acceleration, since the individual results coming from each galaxy are not compatible with a single acceleration. Our findings imply that MOND, although useful for disclosing certain general properties of galaxies, cannot be a fundamental theory, even at the level of galaxies. Consequently, any proposed theory that extends MOND while preserving its dynamics at the level of galaxies cannot be correct (it is rejected at more than 10). Due to its emergent nature, the RAR may be the result of a complex interplay between dark matter and baryons. Particular realizations of this picture are currently being investigated[18, 19, 35].
Methods
0.1 Data.
0.2 Data quality and assumptions in the data.
The data samples that we use are commonly cited and used as the best samples for disc galaxy analyses. Nonetheless, we clarify that they are not raw observational data, since they depend, in particular, on the galaxy surface brightness decomposition into gas, stellar disc and stellar bulge, interpretations on the inclinations, distances, assumptions on axial symmetry and on choices on data binning. When possible we reduce the dependence of our results on such hypotheses; this is done by using nuisance parameters with flat priors, considering different sets of quality cuts and considering two different samples of data (SPARC and THINGS) which use different approaches. Finally, uncontrolled systematics are highly unlikely to appreciably change the strong significance of our results.
0.3 New numerical packages.
To analise the data of all the 175 galaxies from SPARC many times with different hypothesis, we developed two new packages to automatize these procedures. Both packages can be used for other samples, besides SPARC, but are optimized for the latter. They are based on the Wolfram Language and are named MAGMA (for Mathematica Automatized Galaxy Mass Analysis) and mBayes (Bayesian analysis with Mathematica). The MAGMA package has three main functions, and all of them automatically analyse all the galaxies with a single human interaction. The functions have the following purposes: ) one for finding the best fit of a model, by minimizing the with given constraints, ) another function for generating plots for the individual galaxies with the derived fitted parameters, and ) a last function for generating a single table with relevant data from all the fits (including results from mBayes). The best fit is searched using a differential evolution algorithm many times, with different numerical options in parallel, in a twostep procedure, such that the second step refines the first step results. The main purpose of the mBayes package is to perform the Bayesian analysis and find the confidence contours for all the galaxies. It is the most demanding computational part, and it also uses parallel computing. It imports the functions generated by MAGMA together with the values of the parameters’ best fits. The best fits are only used as an initial input for the exploration of the parameter space. The mBayes package generates triangular plots for each galaxy (see Supplementary Materials for two examples) and a table with the numerical values of the 1to5 confidence intervals (which is what is used to generate Figure 1 and Table 1).
0.4 Newtonian acceleration and circular velocity.
The Newtonian acceleration is by definition . The Newtonian potential is defined from the Poisson equation , where is the gravitational constant and the matter density profile. Galaxy matter densities that can be inferred from observations are here decomposed, in accordance with the SPARC and THINGS conventions, into stellar and gaseous components. The stellar component, when necessary, is further decomposed into bulge and disk. Consequently, the total mass profile of a galaxy is decomposed into bulge, disk and gas according to . These are the baryonic contributions, here we consider no dark matter. Due to the linearity of the Poisson equation, for each of these componentes one can derive a corresponding acceleration contribution. Assuming that a disk galaxy is axisymmetric and rotationally supported, the total acceleration experienced by a small radial interval in a galaxy must be oriented towards the galaxy center. Since matter is not spherically distributed in a disk galaxy, individual contributions to the total Newtonian acceleration may be oriented to the opposite direction. These negative contributions are relevant for some galaxies, and are considered in our codes. If the acceleration of a given galaxy component is towards the galaxy center, then for that component, where is the cylindrical radial coordinate. In the context of galaxy rotation curves, it is common to adopt the following convention valid for any direction of the acceleration: , where is the radial component of the acceleration. Hence, galaxy components that reduce the total acceleration at a given radius are represented by a negative velocity at that radius. See the Supplementary Material for examples of galaxy rotation curves.
Both the SPARC and the THINGS databases provide the circular velocities of each baryonic component, and with the sign convention above. The square of the total Newtonian circular velocity is written as
(2) 
Due to the uncertainties on the stellar masstolight ratios, the dimensionless masstolight ratio constants are inserted in the expression above. The provided data on the circular velocity components are for (this value is adopted for the data presentation in order to easy the conversion to any other masstolight ratios).
0.5 Models: MOND with three different interpolating functions.
All the analyses performed here are done considering MOND, but with three different interpolating functions , see equation (1). For the model analysis, it is useful to invert equation (1) and express the acceleration as a function of . It is in this form that the interpolating functions are implemented in the code. For the Standard, Simple and RARinspired interpolating functions, the accelerations read:
(3) 
0.6 Fitted parameters and function.
For all galaxies we consider a tolerance of for the galaxy distance and a tolerance of a factor two on the masstolight ratios (both for the disk and the bulge, when present). For the SPARC sample, the expected values for and are respectively 0.5 and 0.7 in the 3.6 band (these are fixed for all the galaxies). These parameters, together with the constant , compose the fitted parameters. The best fit relative to a given galaxy is defined as the value that maximizes the likelihood function or, equivalently, minimizes the function,
(4) 
where is the number of data points, is the radius at which the velocity was measured with the error , and is the theoretical velocity from MOND (with a given interpolating function). The relation between the theoretical physical acceleration according to MOND and the velocity is .
Details of the Bayesian analysis are in the Supplementary Material.
0.7 Quality cuts.
We perform the following four quality cuts; see the Supplementary Material for the names of the galaxies that were excluded by each quality cut. First, we remove galaxies that have a very high minimum : in these cases MOND is already a bad model – perhaps because of specific observational or dynamical issues – and it does not make sense to include them in the present analysis as the null hypothesis assumes that MOND does work for individual galaxies. A consistent way to achieve the latter is to set a fixed value below which we reject the (null) hypothesis that MOND provides a good fit. Specifically, in the main analysis we decide not to analyze galaxies for which MOND is rejected at the 5 confidence level, that is, we exclude galaxies for which:
(5) 
where is the cumulative distribution with degrees of freedom, is the observed value and is the error function. In this case it is where is the number of data points and the number of fitted parameters, which is 3 or 4 depending on the presence of the bulge. The reference value above is . This selection eliminates 37 galaxies. In order to increase the robustness of our analysis, we also consider the stricter selection criterion of not analyzing galaxies for which MOND is rejected at the 3 confidence level (reference value of ). This eliminates 62 galaxies.
Secondly, we do not include in the analysis of Figure 1 and Table 1 galaxies whose posterior on does not go to zero for so that the regions are ill defined. These galaxies accept very low values of and would not improve the chances of a fundamental value. This second selection eliminates 27 galaxies.
Thirdly, we do not include the 12 SPARC galaxies that are marked with the quality flag (poor quality, as defined by SPARC). Finally, other 10 SPARC galaxies whose inclination is less than are not considered. These last two criteria give the “RAR flag” (), which is 1 for the 153 SPARC galaxies that were considered when evaluating the RAR[1], and zero otherwise. Due to the elimination of the (close to) faceon galaxies, we do not consider galaxy inclination errors since they have minor impact for[11] . The uncertainties on stellar masstolight ratios and galaxy distance are relevant for the dynamical analysis, and these are taken into consideration through the flat priors on and , as described above.
After these four quality cuts, we are left with 100 galaxies. While 79 galaxies are left in the case of the stricter 3 criterion on the minimum .
0.8 Global best fit.
For each galaxy we compute the mean and the variance of the posterior according to:
(6)  
(7) 
where is the (marginalized) posterior on for the galaxy , with where is the number of galaxies that passed the quality cuts discussed above. Then, in order to find the global bestfit value of we introduce the following statistic:
(8) 
which follows a distribution with degrees of freedom, under the approximation that is distributed according to a Gaussian distribution with mean and variance . The minimization of yields the “global best fit” .
0.9 Confidence level in rejecting MOND.
As the degrees of freedom of the defined in Equation (8) is high (), we can approximate the distribution with a Gaussian with mean and standard deviation . We can then estimate the number of at which MOND is rejected according to:
(9) 
Using any of the three interpolating functions for the SPARC sample we find that (main analysis, 5 criterion for the exclusion of galaxies with high ) and (3 criterion).
0.10 Correlations with galactic parameters and observational data robustness.
Observational data on galaxy rotation curves may in principle be more trustworthy for certain values of the galaxy parameters than for others (e.g., some systematic issues with the data may correlate with luminosity, surface brightness, asymptotic circular velocity, …). The quality cuts previously described do not explicitly address issues of this type. In particular, in high luminosity galaxies it is more frequent to find local features in the rotation curves, and these may be either interpreted as necessary features to be reproduced by galaxy models (like the Renzo rule), or as features that appear due to poor modeling of the baryonic matter[13, 11] (e.g., a significative violation of axial symmetry in the observational data, or an error in addressing a variation of the inclination). In the latter case, a local disagreement between the model and the rotation curve data would not be a problem.
In the Supplementary Material we present our analysis of this topic in more detail. We find at most weak correlations that have no significant impact on our results. In particular, by imposing additional quality cuts, such that galaxies of sufficiently high luminosity are not considered, the tension between the data and the existence of a fundamental acceleration scale can be reduced, but the reduction does not change our main results. Even with these additional quality cuts, a fundamental (universal) acceleration scale is still discarded at more than 10.

We thank Stacy McGaugh for clarifications regarding the SPARC sample and comments on a previous version of this paper. This work made use of SPARC “Spitzer Photometry & Accurate Rotation Curves” and of THINGS “The HI Nearby Galaxy Survey”. DCR and VM thank CNPq and FAPES for partial financial support. ADP was supported by the Chinese Academy of Sciences and by the President’s international fellowship initiative, grant no. 2017 VMA0044. ZD thanks the ministry of science, research and technology of Iran for financial support.

DCR and ADP proposed the study. DCR developed the MAGMA package, performed the minimization analysis and contributed to interpretation and design. VM developed the mBayes package, performed the Bayesian analysis and contributed to interpretation and design. ZD carried out the THINGS sample analysis and raised issues that were essential for the beginning of this project. The first draft was written by DCR and VM, and all the authors contributed to its development.

The authors declare that they have no competing financial interests.

Correspondence and requests for materials should be addressed to DCR or VM (emails: davi.rodrigues@cosmoufes.org, marra@cosmoufes.org).
References
 1. McGaugh, S., Lelli, F. & Schombert, J. Radial Acceleration Relation in Rotationally Supported Galaxies. Phys. Rev. Lett. 117, 201101 (2016). 1609.05917.
 2. Lelli, F., McGaugh, S. S., Schombert, J. M. & Pawlowski, M. S. One Law to Rule Them All: The Radial Acceleration Relation of Galaxies. Astrophys. J. 836, 152 (2017). 1610.08981.
 3. Smolin, L. MOND as a regime of quantum gravity. Phys. Rev. D96, 083523 (2017). 1704.00780.
 4. Akerib, D. S. et al. Results from a search for dark matter in the complete LUX exposure. Phys. Rev. Lett. 118, 021303 (2017). 1608.07648.
 5. Aprile, E. et al. First Dark Matter Search Results from the XENON1T Experiment. Phys. Rev. Lett. 119, 181301 (2017). 1705.06655.
 6. Mo, H., van den Bosch, F. & White, S. Galaxy Formation and Evolution (Cambridge University Press, 2010).
 7. Del Popolo, A. Nonbaryonic Dark Matter in Cosmology. Int. J. Mod. Phys. D23, 1430005 (2014). 1305.0456.
 8. Del Popolo, A. & Le Delliou, M. Small scale problems of the CDM model: a short review. Galaxies 5, 17 (2017). 1606.07790.
 9. Liu, J., Chen, X. & Ji, X. Current status of direct dark matter detection experiments. Nature Phys. 13, 212–216 (2017). 1709.00688.
 10. Rodrigues, D. C., del Popolo, A., Marra, V. & de Oliveira, P. L. C. Evidences against cuspy dark matter halos in large galaxies. MNRAS 470, 2410–2426 (2017). 1701.02698.
 11. Lelli, F., McGaugh, S. S. & Schombert, J. M. SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves. AJ 152, 157 (2016). 1606.09251.
 12. Walter, F. et al. THINGS: The H I Nearby Galaxy Survey. AJ 136, 2563–2647 (2008). 0810.2125.
 13. de Blok, W. J. G. et al. HighResolution Rotation Curves and Galaxy Mass Models from THINGS. AJ 136, 2648–2719 (2008). 0810.2100.
 14. Milgrom, M. A modification of the Newtonian dynamics  Implications for galaxies. ApJ 270, 371–389 (1983).
 15. Sanders, R. H. & McGaugh, S. S. Modified Newtonian Dynamics as an Alternative to Dark Matter. Ann. Rev. Astron. Astrophys. 40, 263–317 (2002). astroph/0204521.
 16. Famaey, B. & McGaugh, S. Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions. Living Rev. Rel. 15, 10 (2012). 1112.3960.
 17. Milgrom, M. Road to MOND: A novel perspective. Phys. Rev. D92, 044014 (2015). 1507.05741.
 18. Ludlow, A. D. et al. MassDiscrepancy Acceleration Relation: A Natural Outcome of Galaxy Formation in Cold Dark Matter Halos. Phys. Rev. Lett. 118, 161103 (2017). 1610.07663.
 19. Navarro, J. F. et al. The origin of the mass discrepancyacceleration relation in CDM. Mon. Not. Roy. Astron. Soc. 471, 1841 (2017). 1612.06329.
 20. de Blok, W. J. G. & McGaugh, S. S. Testing Modified Newtonian Dynamics with Low Surface Brightness Galaxies: Rotation Curve FITS. ApJ 508, 132–140 (1998). astroph/9805120.
 21. Gentile, G., Famaey, B., Zhao, H. & Salucci, P. Universality of galactic surface densities within one dark halo scalelength. Nature 461, 627 (2009). 0909.5203.
 22. Tully, R. & Fisher, J. A New method of determining distances to galaxies. AA 54, 661–673 (1977).
 23. McGaugh, S. S., Schombert, J. M., Bothun, G. D. & de Blok, W. The Baryonic TullyFisher relation. ApJ 533, L99–L102 (2000). astroph/0003001.
 24. Hees, A., Famaey, B., Angus, G. W. & Gentile, G. Combined Solar System and rotation curve constraints on MOND. Mon. Not. Roy. Astron. Soc. 455, 449–461 (2016). 1510.01369.
 25. Famaey, B. & Binney, J. Modified Newtonian Dynamics in the Milky Way. MNRAS 363, 603–608 (2005). astroph/0506723.
 26. Gentile, G., Famaey, B. & de Blok, W. THINGS about MOND. AA 527, A76 (2011). 1011.4148.
 27. Dodelson, S. The Real Problem with MOND. Int. J. Mod. Phys. D20, 2749–2753 (2011). 1112.1320.
 28. Sanders, R. H. Clusters of galaxies with modified Newtonian dynamics (MOND). Mon. Not. Roy. Astron. Soc. 342, 901 (2003). astroph/0212293.
 29. Angus, G. W. Are sterile neutrinos consistent with clusters, the CMB and MOND? Mon. Not. Roy. Astron. Soc. 394, 527 (2009). 0805.4014.
 30. Milgrom, M. Bimetric MOND gravity. Phys. Rev. D80, 123536 (2009). 0912.0790.
 31. Babichev, E., Deffayet, C. & EspositoFarese, G. Improving relativistic MOND with Galileon kmouflage. Phys. Rev. D84, 061502 (2011). 1106.2538.
 32. Verlinde, E. P. Emergent Gravity and the Dark Universe. SciPost Phys. 2, 016 (2017). 1611.02269.
 33. Iocco, F., Pato, M. & Bertone, G. Testing modified Newtonian dynamics in the Milky Way. Phys. Rev. D92, 084046 (2015). 1505.05181.
 34. Randriamampandry, T. & Carignan, C. Galaxy mass models: MOND versus dark matter haloes. Mon. Not. Roy. Astron. Soc. 439, 2132–2145 (2014). 1401.5619.
 35. Famaey, B., Khoury, J. & Penco, R. Emergence of the mass discrepancyacceleration relation from dark matterbaryon interactions. JCAP 1803, 038 (2018). 1712.01316.