# On the incompatibility of the Radial Acceleration Relation and Solar System quadrupole in modified gravity MOND

Harry Desmond, Aurélien Hees, Benoit Famaey

► To cite this version:

Harry Desmond, Aurélien Hees, Benoit Famaey. On the incompatibility of the Radial Acceleration Relation and Solar System quadrupole in modified gravity MOND. Monthly Notices of the Royal Astronomical Society, 2024, 530 (2), pp.1781-1795. 10.1093/mnras/stae955. hal-04416917

HAL Id: hal-04416917

https://hal.science/hal-04416917v1

Submitted on 26 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# On the tension between the radial acceleration relation and Solar system quadrupole in modified gravity MOND

Harry Desmond<sup>1,1*</sup> Aurélien Hees<sup>2,2*</sup> and Benoit Famaey<sup>3,3*</sup>

Accepted 2024 April 3. Received 2024 March 26; in original form 2024 January 10

## ABSTRACT

Modified Newtonian dynamics (MOND), postulating a breakdown of Newtonian mechanics at low accelerations, has considerable success at explaining galaxy kinematics. However, the quadrupole of the gravitational field of the Solar system (SS) provides a strong constraint on the way in which Newtonian gravity can be modified. In this paper, we assess the extent to which the AQUAdratic Lagrangian (AQUAL) and QUasilinear MOND (QUMOND) modified gravity formulations of MOND are capable of accounting simultaneously for the radial acceleration relation (RAR), the Cassini measurement of the SS quadrupole and the kinematics of wide binaries in the Solar neighbourhood. We achieve this by inferring the location and sharpness of the MOND transition from the Spitzer Photometry and Accurate Rotation Curves (SPARC) RAR under broad assumptions for the behaviour of the interpolating function and external field effect. We constrain the same quantities from the SS quadrupole, finding that this requires a significantly sharper transition between the deep-MOND and Newtonian regimes than is allowed by the RAR (an $8.7\sigma$ tension under fiducial model assumptions). This may be relieved somewhat by allowing additional freedom in galaxies' mass-to-light ratios – which also improves the RAR fit – and more significantly (to $1.9\sigma$) by removing galaxies with bulges. For the first time, we also apply to the SPARC RAR fit an AQUAL correction for flattened systems, obtaining similar results. Finally, we show that the SS quadrupole constraint implies, to high precision, no deviation from Newtonian gravity in nearby wide binaries, and speculate on possible resolutions of this tension between SS and galaxy data within the MOND paradigm.

**Key words:** gravitation – ephemerides – planets and satellites: general – galaxies: kinematics and dynamics – galaxies: statistics – dark matter.

## 1 INTRODUCTION

The goal of the theory of gravity is to explain as many gravitational phenomena as possible with as few parameters. The theory of general relativity (GR), based on Newtonian gravity, is successful on scales from sub-mm (from Eötvös type experiments; Adelberger, Heckel & Nelson 2003) to tens of Gpc (the dynamics of the cosmos as a whole, assuming dark matter and dark energy). On galaxy scales, the principal alternative is Modified Newtonian dynamics (MOND; Milgrom 1983a,b,c), which postulates a breakdown of Newtonian gravity or inertia at accelerations $\lesssim 10^{-10} \text{ m s}^{-2}$. This has success at explaining several aspects of galaxy dynamics as well as a smattering of observations at other scales (see Famaey & McGaugh 2012; Banik & Zhao 2022, for reviews).

MOND provides a framework within which to understand the otherwise surprising simplicity of galaxy dynamics. For example, the baryonic Tully–Fisher relation – the correlation between the asymptotic circular velocity and baryonic mass of rotation-supported galaxies – has very small intrinsic scatter and an almost perfect power-law shape across five orders of magnitude in mass (McGaugh

et al. 2000; Desmond 2017b). There are no residual correlations with galaxy size or baryonic surface density (Famaey & McGaugh 2012; Lelli, McGaugh & Schombert 2016b; Desmond et al. 2019), which one would typically expect to be produced by the dark matter halos of a Lambda cold dark matter ($\Lambda$CDM) cosmology (based on GR; Desmond & Wechsler 2015). At the same time, galaxies of fixed baryonic mass, while sharing a similar asymptotic circular velocity, display rotation curves (RCs) with a vast diversity of inner shapes even in dwarf galaxies where, in $\Lambda$CDM, dark matter supplies most of the gravity. Indeed, it has long been known that the shapes of RCs *do* depend strongly on baryonic surface density, even in dark matter-dominated galaxies (e.g. Zwaan et al. 1995; de Blok & McGaugh 1997; Swaters et al. 2009). The RC shapes are thus not only diverse at a given mass-scale while sharing a single asymptotic velocity, but are uniform at a given baryonic surface density scale. In $\Lambda$CDM, this implies a diversity of cored and cuspy dark matter profiles which is not currently understood (Oman et al. 2015; Ghari et al. 2019b). This may be explained by highly non-circular motions (Roper et al. 2023), at the price of producing distorted 2D velocity fields not resembling observed ones (e.g. Kuzio de Naray & Kaufmann 2011).

All of this phenomenology for disc galaxies can actually be deduced from a single, apparently more fundamental, local relation (McGaugh, Lelli & Schombert 2016; Lelli et al. 2017) between the

<footnote n="1"><sup>1</sup>Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth PO1 3FX, UK</footnote>

<footnote n="2"><sup>2</sup>SYRTE, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, LNE, 61 avenue de l'Observatoire, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, F-67000 Strasbourg, France</footnote>

<footnote>*E-mail: harry.desmond@physics.ox.ac.uk (HD); aurelien.hees@obspm.fr (AH); benoit.famaey@astro.unistra.fr (BF)</footnote>
total gravitational acceleration at any point in the disc ($g_{\text{obs}}$) and that generated by the observed baryons ($g_{\text{bar}}$). This radial acceleration relation (RAR) has recently been shown to pass all posited tests for fundamentality, and appears to be the root cause of all other radial dynamical correlations of disc galaxies (Stiskalek & Desmond 2023). For such systems, the RAR is nothing more than MOND written in terms of observables, so that, within that paradigm, all the above observations become readily understandable. There were already hints at the inception of MOND that the characteristic rotation velocities of galaxies were correlated with their luminosity (Tully & Fisher 1977); Milgrom (1983a, 1983b, 1983c) promoted them to a tentative law of Nature by proposing that Newtonian dynamics, the weak-field limit of GR, breaks down at ultralow accelerations typical of the outskirts of galaxies (below a new universal constant $a_0 \simeq 10^{-10} \text{ m s}^{-2}$) where one instead has $g_{\text{obs}} = \sqrt{a_0 g_{\text{bar}}}$. In semi-analytic and hydrodynamical models of galaxy formation of $\Lambda$CDM, on the other hand, while the global shape of the RAR can be reproduced, its low scatter and lack of residual correlations are very difficult to understand (e.g. Di Cintio & Lelli 2016; Desmond 2017a; Paranjape & Sheth 2021).

Given MOND's success at the scale of galaxies, it is crucial to investigate the extent to which it may be extended to other regimes. On larger scales, it has problems accounting for the mass discrepancy in clusters (e.g. The & White 1988; Sanders 1999; Pointecouteau & Silk 2005; Angus, Famaey & Buote 2008; Li et al. 2023) and a fully fledged and observationally consistent cosmology is proving challenging to formulate (although see Skordis & Złośnik 2021 for recent progress). On smaller scales, within the Milky Way, key tests are the orbits of wide binary (WB) systems in which the internal orbital acceleration lies below $a_0$, and of bodies in the SS. Both of these systems are embedded within the gravitational field of the Milky Way at or near the position of the Sun, which is around $1.8 a_0$ (Gaia Collaboration et al. 2021). This makes the tests less clean than the fully weak-field probes in the outskirts of low surface density galaxies, but it is nevertheless clear – at least within canonical modified gravity formulations – that MOND should provide a small boost to the dynamics in these systems if it is to provide one for RCs at accelerations of $1.8 a_0$. WBS have had a long and confused history in the context of MOND, with some authors finding consistency with Newtonian dynamics (Pittordis & Sutherland 2023; Banik et al. 2024), while others argue that these analyses are flawed and the data instead prefer MOND (Hernandez et al. 2019, 2023; Hernandez, Cookson & Cortés 2022; Hernandez 2023; Chae 2023a, b; Hernandez & Chae 2023).

In this work, we focus instead on the SS tests, which were revolutionized by the Cassini spacecraft's radioscience tracking (Wolf & Smith 1995; Antreasian et al. 2005; Jacobson et al. 2006; Iess et al. 2010, 2012). Cassini's measurements made it possible to reconstruct the spacecraft's trajectory as it orbited Saturn between 2004 and 2013, and to improve significantly our knowledge of Saturn's orbit. The most constraining measurement in the MOND context is a null detection of the quadrupole moment of the gravitational field of the Sun (Hees et al. 2014). A quadrupole is expected in (at least the most popular) modified gravity formulations of MOND – the AQUAdratic Lagrangian (AQUAL; Bekenstein & Milgrom 1984) and QUasilinear MOND (QUMOND; Milgrom 2010) formulations – given the external field effect (EFE) due to the relative magnitude and orientation of the (internal) Saturn–Sun and (external) Sun–Galactic Centre gravitational fields (Milgrom 2009; Blanchet & Novak 2011). Hees et al. (2016) found Cassini's quadrupole measurement to be inconsistent with a large variety of Newton-to-deep-MOND transitions that describe the RAR well. Although the Cassini mission

ended in 2017, studies of RCs have greatly advanced since then, making the topic ripe for revisiting.

The EFE arises because, as an acceleration-based modification to Newtonian mechanics, MOND is sensitive to the total gravitational field and hence violates the strong equivalence principle (see e.g. Famaey & McGaugh 2012; Banik & Zhao 2022, for a detailed account). While the majority of evidence seems to point towards the helpfulness of the EFE for explaining the kinematics and tidal stability of various types of galaxies within a MOND context (e.g. Famaey, Bruneton & Zhao 2007b; Wu et al. 2008; Haghi et al. 2009, 2016, 2019; McGaugh & Milgrom 2013; Hees et al. 2016; Famaey, McGaugh & Milgrom 2018; Kroupa et al. 2018; Thomas et al. 2018; Banik et al. 2020; Chae et al. 2020; Oria et al. 2021; Asencio et al. 2022; Kroupa et al. 2022), there are also cases where the EFE should be present but appears absent (Freundlich et al. 2022). It is also important to note that, while some form of EFE is generically expected by MOND, in modified inertia formulations it may be very different (e.g. dependent on an object's entire past trajectory) or effectively absent (e.g. Milgrom 2011, 2023a).

The consistency – or lack thereof – between the RAR and SS quadrupole in MOND hinges on the nature of the transition between the Newtonian and deep-MOND regimes, $g_{\text{obs}} = f(g_{\text{bar}})$ in the language of the RAR. The functional form cannot, at present, be derived theoretically and thus must be constrained empirically: the only requirements in MOND (absent the EFE) are the limits $f(g_{\text{bar}}) \to g_{\text{bar}}$ as $g_{\text{bar}} \to \infty$ and $f(g_{\text{bar}}) \to \sqrt{a_0 g_{\text{bar}}}$ as $g_{\text{bar}} \to 0$. Any $f$ that satisfies these limits may be re-expressed as a MONDian 'interpolating function' (IF). Such transition functions are not part of MOND's basic tenets but are used, in one way or another, by all known modified gravity formulations of the theory, generally put by hand into the Lagrangian. The IF could, however, emerge from a fundamental underlying theory with an actually different functional form in different systems, as may in fact be expected in modified inertia formulations (Milgrom 2023a).

Many functions satisfying the MOND requirements have been proposed and investigated in the literature (see Famaey & McGaugh 2012; Banik & Zhao 2022, and references therein). Desmond, Bartlett & Ferreira (2023) evaluated all possible functions of low complexity on dynamical galaxy data, concluding that the optimal function for RC data is fairly complex – and may not have MOND limits – but that more data is required to deduce it unambiguously. In particular, it is not clear from RC data alone that $g_{\text{obs}} \propto \sqrt{g_{\text{bar}}}$ at low $g_{\text{bar}}$. Additional uncertainty arises from the EFE, which depends on the imperfectly observationally characterized baryonic large-scale structure of the Universe and has a functional form that cannot be deduced analytically for general mass distributions. Here, we group IFs into classes or parametric 'families', which contain not only $a_0$ as a degree of freedom but also a *shape* parameter controlling the sharpness of the transition from the deep-MOND to Newtonian regimes.

The aim of this paper is twofold. First, we extend the RAR analysis of Desmond (2023) to cover three IF families under four different state-of-the-art models for the EFE. Within classical modified gravity formulations (AQUAL and QUMOND), this essentially fully spans the space of possibilities both for the form of the MOND force law and the effect of the environment on the galaxies in question, thus affording a definitive inference of the associated parameters with the galaxies' systematically uncertain properties fully accounted for. We also explore various priors for the mass-to-light ratios (M/Ls) of the galaxies' components, in case the fiducial model is overly restrictive, and we explicitly check the effect of the flat disc geometry in the context of AQUAL. Secondly, we compare the RAR {$a_0$,
shape} constraints with those inferred from the measurement of the SS quadrupole by the Cassini mission, and to results of the Banik et al. (2024) wide binary test (WBT). These all provide independent probes of MOND, which, if inconsistent with each other, would pose a severe problem for the AQUAL (Bekenstein & Milgrom 1984) and QUMOND (Milgrom 2010) weak-field modified gravity formulations.

In Section 2, we describe the galaxy and SS data that we use as constraints. In Section 3, we lay out our MOND models (Section 3.1), our method for inferring their parameters from the RAR (Section 3.2) and SS quadrupole (Section 3.3), and the method for relating these to the WBT (Section 3.4). Section 4 describes the results, separately for the RAR (Section 4.1) and SS quadrupole (Section 4.2), and explores why the two sets of constraints are difficult to reconcile (Section 4.3). We conclude in Section 5. Throughout the paper, log has base 10.

# 2 OBSERVATIONAL DATA

## 2.1 SPARC galaxy sample

For the RAR, we utilize the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample (Lelli, McGaugh & Schombert 2016a),<fnref n="1" /> containing 175 RCs from the literature with photometry at 3.6 µm from the Spitzer satellite (see also Gentile, Famaey & de Blok 2011). We adopt the same quality cuts as Lelli et al. (2017), excluding galaxies with quality flag 3 (indicating strong asymmetries, non-circular motions and/or offsets between the stellar and H1 distributions) or inclination $i < 30^{\circ}$, and points with a fractional rotation velocity uncertainty >10 per cent. 2696 points from 147 galaxies remain, of which 31 contain a central stellar bulge. This is the same sample as was used in Desmond (2023).

## 2.2 Cassini measurement

The MOND paradigm breaks the strong equivalence principle such that the external gravitational field of our Galaxy impacts the internal dynamics of the SS (Milgrom 2009; Blanchet & Novak 2011). The leading effect is a modification to the central Newtonian potential of the Sun by an additional quadrupole term (adopting the Einstein summation convention)

$$ \delta\Phi(\mathbf{x}) = -\frac{Q_2}{2} x^i x^j \left( \hat{e}_i \hat{e}_j - \frac{1}{3} \delta_{ij} \right), \qquad (1) $$

where $\hat{e} = g_{ext}/g_{ext}$ is a unit vector pointing towards the Galactic Centre, $x$ the position within the SS with respect to the Sun, $\delta_{ij}$ the Kronecker delta and $Q_2$ a parameter that depends on the MOND IF and acceleration scale $a_0$, the Sun's mass $M$ and the value of the external gravitational field from the Galaxy, $g_{ext}$.

The modified Newtonian potential from equation (1) will induce an anomalous acceleration which rises linearly with distance from the Sun and can be sought using the kinematics of planets (Milgrom 2009; Blanchet & Novak 2011; Hees et al. 2012), Kuiper Belt objects (Brown & Mathur 2023), asteroids, or comets (Maquet & Pierret 2015; Vokrouhlický, Nesvorný & Tremaine 2024). In Hees et al. (2014), the $Q_2$ parameter of equation (1) was inferred using the DE430 planetary ephemerides data (Folkner et al. 2014). Of central importance for constraining this parameter are the 9 yr of Cassini range and Doppler tracking data which strongly constrain Saturn's

orbit. Although subject to some systematic uncertainty (see Hees et al. 2014 for details), the $1\sigma$ constraint may be described by

$$ Q_2 = (3 \pm 3) \times 10^{-27} \text{ s}^{-2}. \qquad (2) $$

# 3 METHOD

## 3.1 Parametrising MOND

We consider a diverse set of possible MONDian descriptions of the relation between total and baryonic acceleration, viz. IFs. The MOND force law in highly symmetric configurations can be expressed as

$$ g = g_N \nu \left( \frac{g_N}{a_0} \right), \qquad (3) $$

where $g$ is the total gravitational field (or dynamical acceleration), $g_N$ is the Newtonian gravitational field sourced by baryons, $\nu(y)$ is the IF and $g_N \equiv |g_N|$. The only a priori requirement on $\nu$ from the basic tenets of MOND is the asymptotic behaviour $\nu(y) \to 1$ when $y \to \infty$, and $\nu(y) \to 1/\sqrt{y}$ when $y \to 0$. The three most common examples are the ‘Simple’ IF (Famaey & Binney 2005; Zhao & Famaey 2006)

$$ \nu_{\text{simp}}(y) = \frac{1 + (1 + 4y^{-1})^{1/2}}{2}, \qquad (4) $$

the ‘Standard’ IF (Milgrom 1983c):

$$ \nu_{\text{stand}}(y) = \left( \frac{1 + (1 + 4y^{-2})^{1/2}}{2} \right)^{1/2} \qquad (5) $$

and the ‘RAR’ (also called ‘McGaugh–Lelli–Schombert’ or ‘MLS’) IF (McGaugh et al. 2016)

$$ \nu_{\text{RAR}}(y) = (1 - \exp(-y^{1/2}))^{-1}. \qquad (6) $$

Absent theoretical or compelling observational evidence for a particular IF, it is convenient to group IFs into ‘families’ sharing a more general functional form in order to test them. These have a parameter (in addition to $g_N/a_0$, their only other degree of freedom) that interpolates between specific IFs by controlling the sharpness of the transition between the deep-MOND ($g_N/a_0 \ll 1$) and Newtonian ($g_N/a_0 \gg 1$) regimes. Here, we investigate three IF families described in Famaey & McGaugh (2012):

$$ \nu_n(x) = \left[ \frac{1 + (1 + 4x^{-n})^{1/2}}{2} \right]^{1/n}, \qquad (7a) $$

$$ \nu_{\delta}(x) = (1 - e^{-x^{\delta/2}})^{-1/\delta}, \qquad (7b) $$

$$ \nu_{\gamma}(x) = (1 - e^{-x^{\gamma/2}})^{-1/\gamma} + (1 - 1/\gamma)e^{-x^{\gamma/2}}. \qquad (7c) $$

These cover the great majority of functions that have been considered in the literature; in particular, the $n$-family encompasses the Simple ($n=1$) and Standard ($n=2$) IFs, while the $\delta$- and $\gamma$-families contain the RAR IF at $\delta = \gamma = 1$.<fnref n="2" /> We refer to $n$, $\delta$ and $\gamma$ collectively as ‘shape’.

As well as a model without any EFE, we explore the possible consequences of the EFE on the RAR fit using two different formulae. These are the ‘Freundlich–Oria analytic’ (Oria et al. 2021; Freundlich et al. 2022) and ‘AQUAL numerical’ models of Chae & Milgrom

<footnote n="2">²There is also a $\beta$ family (equation 51 of Famaey & McGaugh 2012), but this is not flexible enough to cover Simple- or RAR-like behaviour, and hence is poor at fitting the RAR.</footnote>

<footnote n="1">¹http://astroweb.cwru.edu/SPARC/</footnote>
(2022) (cases 4 and 6 in their table 1). These are both fitting functions for RCs, the former developed in the context of the QUMOND model and the latter that of AQUAL. These modify the raw IFs and hence can be applied to any of them. Considering the Newtonian external field strength in units of $a_0$, i.e. $e_N \equiv g_{N, ext}/a_0$, the QUMOND EFE model is given by

$$v_{\text{EFE, QUMOND}}(y) = v \left( \min \left[ y + \frac{e_N^2}{3y}, e_N + \frac{y^2}{3e_N} \right] \right) \quad (8)$$

and the AQUAL model by

$$v_{\text{EFE, AQUAL}}(y) = v(y_{\beta}) \left[ 1 + \tanh \left( \frac{\beta e_N}{y} \right)^{\gamma} \frac{\hat{v}(y_{\beta})}{3} \right], \quad (9)$$

where $\hat{v}(y) \equiv d \ln v(y)/d \ln y$ and $y_{\beta} \equiv \sqrt{y^2 + (\beta e_N)^2}$. The best-fitting values of the parameters $\beta$ and $\gamma$ are 1.1 and 1.2, respectively (Chae & Milgrom 2022). These roughly span the space of possible EFE behaviour, although note that equation (9) is strictly only intended to be used in the outer parts of RCs (where the EFE plays a larger role).

For each of these models, we then consider two different ways of setting the external field strengths $e_N$ of the SPARC galaxies. The first is simply to treat $e_N$ as a global constant and infer it (with a wide uniform prior) from the data. The second is to allow the local $e_N$ to vary galaxy-by-galaxy, with a prior specified by the ‘maximum clustering’ model of Chae et al. (2021), which infers $e_N$ for each SPARC galaxy from its large-scale baryonic environment assuming unseen baryons to correlate maximally with observed baryons.<fnref n="3" /> $e_N$ for each galaxy separately is then a free parameter of the inference to be marginalized over. As in Desmond (2023), for galaxies outside the footprint of the Sloan Digital Sky Survey (where Chae et al. 2021 could not calculate precise $e_N$ values), we use the median over all SPARC galaxies within the footprint as the prior centre, with an uncertainty twice the median uncertainty for those galaxies. For the maximum-clustering case, this gives the prior $\log(e_N) = -2.300 \pm 0.575$.

Finally, it is important to note that in MOND equation (3) (often referred to as the ‘algebraic MOND relation’) is exact only for circular orbits in modified inertia versions of the paradigm (Milgrom 1994, 2022), and cannot be strictly exact in modified gravity, even in the absence of an EFE. We will consider it a sufficient approximation (Jones-Smith et al. 2018; Oria et al. 2021) in our fiducial analysis, but in Section 4.1.3, we will investigate an alternative formula for the MOND boost taking into account the correction for flattened systems within the AQUAL formulation.

## 3.2 RAR inference

Our method to analyse the RAR extends that of Desmond (2023, see especially sec. 3.2). In short, we infer the parameters of the RAR ($a_0$ and intrinsic scatter $\sigma_{\text{int}}$, in addition to the IF shape) simultaneously with the parameters describing each galaxy using priors from previous measurements and theoretical expectations. The galaxy nuisance parameters are distance $D$, inclination $i$,

luminosity $L_{3,6}$, the M/L of the disc $\Upsilon_{\text{disc}}$, bulge $\Upsilon_{\text{bulge}}$ and gas $\Upsilon_{\text{gas}}$, and the external field strengths $e_N$. This amounts to ~900 parameters in total, which is too high a dimensionality for most Markov chain Monte Carlo samplers to handle reliably. It is, however, necessary in order to properly propagate the galaxies’ parameters as systematic rather than statistical uncertainties, and map out their degeneracies with the global properties of the RAR. Leveraging automatic differentiation in Jax, we employ the No U-Turns sampler (a species of Hamiltonian Monte Carlo) as implemented in `numpysro` (Bingham et al. 2019; Phan, Pradhan & Jankowiak 2019). This produces a fully converged chain of ~2000 points in ~10 min; we concatenate 28 chains with multiprocessing for improved statistics and to verify insensitivity to the initialization. Restricting to the Simple IF and a single parametrization of the EFE, Desmond (2023) revealed the intrinsic scatter of the RAR to be minute ($\sigma_{\text{int}} = 0.034 \pm 0.001(\text{stat}) \pm 0.001(\text{sys})$) dex, supporting the claim that the RAR is ‘tantamount to a law of nature’ (Lelli et al. 2017). Weak evidence was adduced for the EFE.

Besides freeing up the IF and EFE implementations as discussed above, we now develop more flexible models for the M/L of the SPARC galaxies’ discs and bulges. These are the only parameters of the galaxies that are not directly constrained empirically, yet are crucial for locating the galaxies’ RC points on the RAR plane and hence constraining the RAR parameters. In particular, we consider six possible priors:

(i) The fiducial SPARC model in which $\Upsilon_{\text{disc}}$ and $\Upsilon_{\text{bulge}}$ follow lognormal<fnref n="4" /> priors with means 0.5 and 0.7 and widths 0.125 and 0.175, respectively (this is also the model used in Desmond 2023).

(ii) $\Upsilon_{\text{disc}}$ and $\Upsilon_{\text{bulge}}$ drawn from separate Gaussian hyperpriors with means $\mu_d$ and $\mu_b$ inferred from the data with independent wide uniform priors. We retain a width of 25 per cent for the hyperprior. This models a scenario in which the centres of the prior M/L distributions are unknown, although their uncertainties follow the SPARC error model.

(iii) As (ii) but requiring $\mu_b > \mu_d$, as expected from population synthesis models.

(iv) $\Upsilon_{\text{disc}}$ and $\Upsilon_{\text{bulge}}$ drawn independently from wide uniform priors, modelling a scenario in which nothing is known about them a priori.

(v) As (iv) but requiring $\Upsilon_{\text{bulge}} > \Upsilon_{\text{disc}}$ on a galaxy-by-galaxy basis, as expected from population synthesis models.

(vi) As (ii) but removing the 31 galaxies with bulges. This models a scenario in which bulges are too poorly understood to be included in the sample.

(vii) As (ii) but removing the 116 galaxies without bulges. This provides a counterpoint to (vi), allowing us to investigate the consistency of bulge and bulge-free galaxies.<fnref n="5" />

The $\Upsilon_{\text{gas}}$ model, based on McGaugh et al. (2020) with a 10 per cent uncertainty, is not altered. Adding extra dark gas has been shown not to improve MOND RC fits (Ghari, Haghi & Zonoozi 2019a).

We use uniform priors on $a_0$, $\sigma_{\text{int}}$, and shape, all of which are well enough constrained by the data for the choice of prior to be unimportant. In particular, using a log-prior on the dimensionful quantities $a_0$ or $\sigma_{\text{int}}$, or a prior on shape corresponding to a flat prior

<footnote n="3">³Chae et al. (2021) also obtain ‘no clustering’ results in which unseen baryons do not correlate at all with observed baryons and hence do not systematically increase $e_N$ but only its scatter, and Desmond (2023) considers an ‘average clustering’ model midway between the two. We focus on the maximum-clustering case because it was found to yield best agreement with the data in Chae et al. (2021, 2022), and is a priori the most likely case in MOND. The differences for our purposes are minor.</footnote>

<footnote n="4">⁴We use lognormal for consistency with previous works, but the results with a normal prior are almost identical.</footnote>

<footnote n="5">⁵Very similar results for models (vi) and (vii) are obtained using free uniform priors rather than Gaussian hyperpriors for the remaining $\Upsilon$ parameters.</footnote>
on $Q_2$ (see Section 3.3) makes negligible difference to the RAR results.

### 3.3 Quadrupole inference

Theoretically, the SS quadrupole value that appears in equation (1) can be expressed as (Milgrom 2009)

$$Q_2 = -\frac{3a_0^{3/2}}{2\sqrt{GM}}q(\tilde{e}), \quad (10)$$

where, in the SS, $M = 1 M_\odot$, and $q$ is a dimensionless parameter that depends only on the MOND IF as well as the value of

$$\tilde{e} \equiv \frac{g_{\text{ext}}}{a_0}, \quad (11)$$

where $g_{\text{ext}}$, in this case, is the external field of the Milky Way at the SS. In the context of QUMOND, Milgrom (2009) has derived an exact expression for $q$:

$$q(\tilde{e}) = \frac{3}{2} \int_{0}^{\infty} dv \int_{-1}^{1} d\xi (v-1)[e_N(3\xi - 5\xi^3) + v^2(1-3\xi^2)], \quad (12)$$

where $v = v[\sqrt{e_N^2 + v^4 + 2e_N v^2 \xi}]$ and $e_N$ is the solution of $e_N v (e_N) = \tilde{e}$. In the case of AQUAL, the above integral leads only to an approximate value for the $q$ parameter. In this case, $q$ must be computed by numerically solving the non-linear Poisson equation (as done in Milgrom 2009 and Blanchet & Novak 2011). By comparing the QUMOND values with those obtained using AQUAL in table 1 of Milgrom (2009) and table 1 of Blanchet & Novak (2011), it is clear that our QUMOND calculation leads to a lower $Q_2$,$^6$ and is therefore conservative when it comes to assessing the tension with the data.

The value of $q(\tilde{e})$ obtained from equation (12) depends mainly on the behaviour of the IF around $e_N$. To illustrate this, Fig. 1 presents, as a function of the argument of the $v$-function in a small bin between $y_0$ and $y_0 + \Delta y$ (with $\Delta y$ set to 0.1), the contribution $\Delta q(y_0)$ of $v(y)$ in this bin to the value of $q(\tilde{e})$. In practice, $\Delta q(y_0)$ is computed by replacing the function $v(y)$ by 1 in equation (12) (i.e. $g = g_N$) except in the range $[y_0, y_0 + \Delta y]$ where it instead follows equation (6) (i.e. $\delta = 1$ in equation 7b or $\gamma = 1$ in equation 7c). Fig. 1 presents the variation of $\Delta q$ with $y_0$ for three different values of the external field $\tilde{e}$. The vertical-dashed lines show the corresponding values of $e_N$. It can be noticed that $|\Delta q|$ is large only around $e_N$, which is the case for all other IFs too. This means that the Cassini constraint probes the IF around $y \sim e_N = g_{N, ext}/a_0$, the gravitational field of the Milky Way at the Sun, similarly to the WBT (Banik et al. 2024) and to the analysis of Vokrouhlický et al. (2024). This is independent of whether the internal Newtonian acceleration probed is higher (as in the Cassini constraint) or lower (as in the aphelia of long-period comets) than $g_{N, ext}$ in these different cases.

Then, for each IF family (equation 7), we have computed the $q$ factor on a regular grid for the ($\tilde{e}$, shape) parameters with $\tilde{e}$ ranging between 1 and 20 (with spacing 0.05) and the shape parameter between 0.5 and 10.4 (with spacing 0.1). To reach lower values of $a_0$, we supplement this with an irregular sampling of $\tilde{e}$ values extending to 1000. We then interpolate this grid using a 2D cubic spline. We use equation (2) as the likelihood function in our inference

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Variation of $\Delta q(y_0)$, the contribution from $v(y)$ between $y_0$ and $y_0 + 0.1$ to $q(\tilde{e})$, for the IF $v_{\text{RAR}}$ of equation (6) and three different values of $\tilde{e}$. The value of $q(\tilde{e})$ from equation (12) equals the sum of the points for each curve, which gives $q(1) = 0.094$, $q(1.5) = 0.159$, and $q(2) = 0.221$. The vertical dashed lines represent the corresponding values of $e_N$. The plot shows that the value of $q(\tilde{e})$ mainly probes the behaviour of the IF around $e_N$.</figcaption></figure>

to constrain $\{a_0, \text{shape}\}$. Our model for the external field of the SS sourced by the Milky Way derives from the Gaia EDR3 measurement of the acceleration of the Sun, which is $g_{\text{ext}} = 2.32 \pm 0.16 \times 10^{-10}$ m s$^{-2}$ (Gaia Collaboration et al. 2021). While the $2\sigma$ lower bound $2 \times 10^{-10}$ m s$^{-2}$ is realistic (corresponding to a large but acceptable peculiar velocity of the Sun; Bovy et al. 2015), the $2\sigma$ upper bound is not (as it would imply an unrealistic negative peculiar velocity of the Sun). We therefore consider the upper bound on $g_{\text{ext}}$ to be $2.48 \times 10^{-10}$ m s$^{-2}$. We treat it in one of two ways: (i) use the value within the range $[2, 2.48] \times 10^{-10}$ m s$^{-2}$ that leads to the lowest predicted value of $Q_2$ in order to be conservative in our inference of $\{a_0, \text{shape}\}$ and maximize consistency with the RAR, or (ii) infer $g_{\text{ext}}$ along with $\{a_0, \text{shape}\}$ using the prior $N(2.32, 0.16)$ truncated at 2.48. These methods give very similar results; we present results for the second method because it is more statistically principled and makes the $Q_2$ and RAR inferences fully independent.

For the $Q_2$ inference, a flat prior on shape produces a posterior extending to infinity, in which limit $Q_2 \to 0$. As well as complicating the construction of confidence intervals, this would introduce a strong volume effect into the PPD of the WBT test statistic (see next subsection). We therefore instead adopt a prior flat on $Q_2$ (with a uniform prior on $a_0$), which matches that used in Hees et al. (2014) to derive the constraint of equation (2). This is achieved by taking the numerical derivative of $Q_2$ with respect to shape at each point within the grid.

### 3.4 Connection to the wide binary test

Following Banik et al. (2024), we summarize the expected result of a WBT in the Solar neighbourhood using a parameter $\alpha_{\text{grav}}$ defined as

$$\alpha_{\text{grav}} \equiv \frac{\sqrt{\eta} - 1}{0.193}, \quad (13)$$

where $\eta \equiv (g_r)/g_N$, the boost of the azimuthally averaged asymptotic radial gravity compared to Newton. In QUMOND, this is given by

$$\eta = v_e \left( 1 + \frac{1}{3} \frac{\partial \ln v_e}{\partial \ln(g_{N, \text{ext}})} \right), \quad (14)$$

<footnote>$^6$More precisely, we compare the values of $Q_2$ obtained for $v_n$ with $n = 1, 2, 5$, and 20. The value of $Q_2$ is always larger in AQUAL, by less than $10^{-27}$ s$^{-2}$, except for $n = 1$ where it is about 25 per cent larger. The tension is therefore if anything underestimated in our QUMOND calculation.</footnote>
where $v_e \equiv v(e_N)$ is the IF evaluated at the Newtonian-equivalent Galactic gravity $e_N = g_{N, ext}/a_0$ at the Solar position. In AQUAL, it is instead given by

$$ \eta = \mu_e^{-1} \left( \frac{\partial \ln \mu_e}{\partial \ln g_{\text{ext}}} \right)^{-1/2} \times \tan^{-1} \left( \left( \frac{\partial \ln \mu_e}{\partial \ln g_{\text{ext}}} \right)^{1/2} \right). \quad (15) $$

$\mu$ is an alternative description of the IF, defined by $\mu(x) = v(y)^{-1}$ where $y = x\mu(x)$. equation (13) is chosen such that $\alpha_{\text{grav}} = 0$ for fully Newtonian gravity, and $\alpha_{\text{grav}} = 1$ for the case of QUMOND with the Simple IF, assuming $a_0 = 1.2 \times 10^{-10} \text{ m s}^{-2}$ and $\tilde{e} = 1.8$. Table 4 of Banik et al. (2024) shows the $\alpha_{\text{grav}}$ values corresponding to a few choices of IF at fixed $a_0$ and $g_{N, ext}$.

Equation (13) describes a mapping from $\{a_0, \text{shape}, g_{N, ext}\}$ to $\alpha_{\text{grav}}$, allowing us to convert results on these parameters from either the RAR or $Q_2$ into a PPD of $\alpha_{\text{grav}}$ values. This is the distribution we would expect to see given the RAR or $Q_2$ inference, and may subsequently be compared to the posterior of a WBT such as that of Banik et al. (2024) to assess the consistency of the MOND interpretation of these systems. For simplicity, following Banik et al. we fix $g_e = 2.142 \times 10^{-10} \text{ m s}^{-2}$ for this comparison, since varying this within the range allowed by Gaia Collaboration et al. (2021) makes little difference to the results. Although straightforward, we do not compute $\{a_0, \text{shape}\}$ posteriors from the inference of Banik et al., as we do in Section 4.2 from the $Q_2$ measurement. This is in order not to restrict ourselves to those results, reflecting the less established nature of the WBT than the $Q_2$ measurement and the fact that other WBT analyses have reached substantially different conclusions (Hernandez et al. 2019; 2022, 2023; Hernandez 2023; Chae 2023a,b).

Like $Q_2$, $\alpha_{\text{grav}}$ goes to 0 in the limit shape $\to \infty$, so a flat prior on shape and the associated semi-infinite prior volume towards that limit would have made $\alpha_{\text{grav}}$ appear more strongly constrained by Cassini than it in fact is. This is cured by adopting a flat prior on $Q_2$. Note that the $\alpha_{\text{grav}}$ inference of Banik et al. (2024), to which we compare the $\alpha_{\text{grav}}$ PPD from $Q_2$, used instead a uniform prior on $\alpha_{\text{grav}}$. We have checked that switching to this prior – or a log-uniform prior in $Q_2$ or $\alpha_{\text{grav}}$ – makes little difference to our results.

# 4 RESULTS

## 4.1 Constraints and predictions from the RAR

### 4.1.1 Fiducial mass-to-light ratio priors

We begin with the fiducial SPARC M/L model (Schombert & McGaugh 2014). Table 1 shows, in this case, the constraints on the RAR parameters for each IF family and EFE model. We do not show $\sigma_{\text{int}}$, which is between 0.033 and 0.035 dex for all models with an uncertainty of 0.001. For the local EFE models, where $e_N$ varies per galaxy, we show the median and 68 per cent confidence intervals of the median $e_N$ values across the posteriors of all galaxies; these entries are italicized to distinguish them from the qualitatively different global $e_N$ constraints appearing in the same column. For each model, we also summarize the PPDs of $Q_2$ and $\alpha_{\text{grav}}$ by their median and 68 per cent confidence interval, including the tension with the Cassini measurement in the case of $Q_2$. This is calculated assuming a Gaussian distribution of the lower uncertainty.

The final two columns describe the goodness-of-fit, specifically the Bayesian information criterion (BIC) relative to the reference RAR IF shown in the top row. The BIC is the limit of the Bayesian evidence (the probability of the data given the model) when the

posterior is modelled as a Gaussian around the maximum a posteriori point and the number of data points greatly exceeds the number of free parameters. Although neither of these assumptions are manifestly true in our case, the BIC still functions as a useful model comparison heuristic by trading off the accuracy of a model with its complexity in terms of number of free parameters. It is given by (Schwarz 1978)

$$ \text{BIC} \equiv k \ln(N) - 2 \ln(\hat{\mathcal{L}}), \quad (16) $$

where $k$ is the number of free parameters, $N$ is the number of data points, and $\hat{\mathcal{L}}$ is the maximum-likelihood value. The evidence is proportional to $\exp(-\text{BIC}/2)$. On the Jeffreys scale (Jeffreys 1939), an evidence ratio in excess of 100 ($|\Delta\text{BIC}| > 9.21$) indicates ‘decisive’ evidence in favour of the higher evidence, lower BIC model.

In the limit of much data, the likelihood approximates the posterior because the prior does not scale with the number of data points while the likelihood does. In our case, however, a better estimator of model quality may be achieved by replacing the maximum-likelihood with maximum-posterior value. This leads us to define

$$ \text{BIC(P)} \equiv k \ln(N) - 2 \ln(\hat{\mathcal{P}}), \quad (17) $$

for maximum posterior value $\hat{\mathcal{P}}$, which we also show. Note that since all log-probability values are negative, introducing additional parameters must increase BIC(P) relative to BIC. For models with the same priors, however, BIC(P) may provide a better model comparison statistic because it is the relative posterior probabilities of the models that are important, not only the likelihoods they assign the data.

We highlight three results from this investigation:

(i) $a_0$ and shape are in all cases well constrained and largely consistent between the models. They show the RAR to have a transition location consistent with recent studies (e.g. table 3 of Desmond 2023) and a transition sharpness closely approximating the Simple (shape = 1 in the n-family) and RAR (shape = 1 in the $\delta$- and $\gamma$-families) IFs. There is therefore at most weak evidence for a more general IF relative to Simple or RAR. A corollary is that the family with which one extends these IFs is not particularly important, allowing us to focus our remaining analyses on a single one. We choose the $\delta$-family because it is the simplest that includes the most popular RAR function. Some historical context to our $a_0$ constraints may be found in Section 5.

(ii) Allowing local gravitational field strengths in the EFE is strongly disfavoured by the BIC relative to the no-EFE model due to the addition of 147 free parameters, and allowing a single global field strength is mildly disfavoured. Our analysis does not therefore show strong evidence for the existence of the EFE in the SPARC RAR. This can also be seen from the fact that the $e_N$ values of the global-EFE models are all consistent with 0 within $2\sigma$. It is, however, curious to note that while the AQUAL constraints on a global $e_N$ are significantly below the expectation from the maximum-clustering prior (and hence also the galaxy-by-galaxy $e_N$ constraints from the local EFE model), the QUMOND constraints are almost identical. This suggests that the RC constraints accord better with the large-scale structure priors in QUMOND than AQUAL, which is supported by the lower BIC and BIC(P) values for the QUMOND EFE models. These results differ from those of Chae et al. (, 2021, 2022) due to their focus on specific galaxies and ours on the overall RAR, and our different fitting and goodness-of-fit assessment procedures.

(iii) $Q_2$ and $\alpha_{\text{grav}}$ are clearly non-zero in all cases, indicating that deviation from Newtonian gravity in the SS quadrupole and WB
<table><caption>**Table 1.** Constraints on RAR (or, equivalently, algebraic MOND) parameters and goodness-of-fit statistics for the fiducial SPARC *M/L* model. The parameter *a*<sub>0</sub> has units of 10<sup>-10</sup> m s<sup>-2</sup> and *Q*<sub>2</sub> of 10<sup>-27</sup> s<sup>-2</sup>. We show all combinations of IF family (equation 7) and EFE model (equations 8 and 9). For the models with galaxy-specific ('local') EFE, the quoted *e*<sub>N</sub> constraints – written in italics – describe the median stacked posteriors over all galaxies, implicitly including the maximum-clustering prior. The *Q*<sub>2</sub> (equation 10) and *α*<sub>grav</sub> (equation 13) columns summarize the posterior predictive distributions (PPDs) of the quadrupole and WB test statistic from the RAR fits, while *σ*<sub>Q2</sub> is the number of sigma tension between the *Q*<sub>2</sub> prediction and the Cassini measurement. For reference, the first row shows the RAR IF fit without the EFE; the BIC values (defined using either the maximum-likelihood or maximum-posterior points) are shown relative to this.</caption><thead><tr><th>IF family</th><th>EFE model</th><th>Shape</th><th>a<sub>0</sub></th><th>e<sub>N</sub></th><th>Q<sub>2</sub></th><th>σ<sub>Q2</sub></th><th>α<sub>grav</sub></th><th>ΔBIC</th><th>ΔBIC(P)</th></tr></thead><tbody><tr><td>RAR IF</td><td>No EFE</td><td>–</td><td>1.03<sup>+0.03</sup><sub>-0.03</sub></td><td>–</td><td>29.2<sup>+0.3</sup><sub>-0.4</sub></td><td>8.7</td><td>0.74<sup>+0.03</sup><sub>-0.03</sub></td><td>0</td><td>0</td></tr><tr><td>δ</td><td>No EFE</td><td>0.97<sup>+0.04</sup><sub>-0.04</sub></td><td>1.02<sup>+0.04</sup><sub>-0.04</sub></td><td>–</td><td>29.4<sup>+0.4</sup><sub>-0.5</sub></td><td>8.7</td><td>0.78<sup>+0.06</sup><sub>-0.06</sub></td><td>11.1</td><td>-10.2</td></tr><tr><td>δ</td><td>AQUAL global</td><td>0.98<sup>+0.04</sup><sub>-0.04</sub></td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0017<sup>+0.001</sup><sub>-0.001</sub></td><td>29.4<sup>+0.4</sup><sub>-0.5</sub></td><td>8.7</td><td>0.77<sup>+0.07</sup><sub>-0.06</sub></td><td>18.7</td><td>15.2</td></tr><tr><td>δ</td><td>AQUAL local</td><td>1.14<sup>+0.05</sup><sub>-0.05</sub></td><td>1.25<sup>+0.05</sup><sub>-0.05</sub></td><td>0.0048<sup>+0.0082</sup><sub>-0.0020</sub></td><td>30.2<sup>+0.5</sup><sub>-0.5</sub></td><td>8.9</td><td>0.71<sup>+0.06</sup><sub>-0.06</sub></td><td>1090</td><td>1480</td></tr><tr><td>δ</td><td>QUMOND global</td><td>0.98<sup>+0.04</sup><sub>-0.04</sub></td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0049<sup>+0.0015</sup><sub>-0.0019</sub></td><td>29.4<sup>+0.4</sup><sub>-0.5</sub></td><td>8.7</td><td>0.77<sup>+0.07</sup><sub>-0.06</sub></td><td>13.9</td><td>-1.98</td></tr><tr><td>δ</td><td>QUMOND local</td><td>1.08<sup>+0.05</sup><sub>-0.05</sub></td><td>1.17<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0050<sup>+0.0030</sup><sub>-0.0018</sub></td><td>29.9<sup>+0.5</sup><sub>-0.5</sub></td><td>8.8</td><td>0.73<sup>+0.07</sup><sub>-0.06</sub></td><td>1090</td><td>1470</td></tr><tr><td>γ</td><td>No EFE</td><td>1.03<sup>+0.07</sup><sub>-0.07</sub></td><td>1.03<sup>+0.03</sup><sub>-0.03</sub></td><td>–</td><td>29.1<sup>+0.4</sup><sub>-0.4</sub></td><td>8.6</td><td>0.71<sup>+0.07</sup><sub>-0.07</sub></td><td>5.93</td><td>-1.99</td></tr><tr><td>γ</td><td>AQUAL global</td><td>1.04<sup>+0.07</sup><sub>-0.07</sub></td><td>1.04<sup>+0.03</sup><sub>-0.03</sub></td><td>0.0018<sup>+0.0010</sup><sub>-0.0010</sub></td><td>29.2<sup>+0.4</sup><sub>-0.4</sub></td><td>8.7</td><td>0.71<sup>+0.07</sup><sub>-0.07</sub></td><td>14.7</td><td>14.9</td></tr><tr><td>γ</td><td>AQUAL local</td><td>1.14<sup>+0.06</sup><sub>-0.06</sub></td><td>1.19<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0048<sup>+0.0070</sup><sub>-0.0020</sub></td><td>30.8<sup>+0.4</sup><sub>-0.4</sub></td><td>9.2</td><td>0.76<sup>+0.06</sup><sub>-0.06</sub></td><td>1100</td><td>1510</td></tr><tr><td>γ</td><td>QUMOND global</td><td>1.04<sup>+0.07</sup><sub>-0.07</sub></td><td>1.04<sup>+0.03</sup><sub>-0.03</sub></td><td>0.0050<sup>+0.015</sup><sub>-0.0018</sub></td><td>29.2<sup>+0.4</sup><sub>-0.4</sub></td><td>8.6</td><td>0.71<sup>+0.07</sup><sub>-0.06</sub></td><td>9.69</td><td>8.32</td></tr><tr><td>γ</td><td>QUMOND local</td><td>1.11<sup>+0.07</sup><sub>-0.07</sub></td><td>1.14<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0050<sup>+0.0030</sup><sub>-0.0017</sub></td><td>30.2<sup>+0.4</sup><sub>-0.4</sub></td><td>9.0</td><td>0.73<sup>+0.06</sup><sub>-0.06</sub></td><td>1080</td><td>1490</td></tr><tr><td>n</td><td>No EFE</td><td>1.02<sup>+0.04</sup><sub>-0.04</sub></td><td>1.08<sup>+0.04</sup><sub>-0.04</sub></td><td>–</td><td>28.4<sup>+0.4</sup><sub>-0.4</sub></td><td>8.4</td><td>0.87<sup>+0.07</sup><sub>-0.06</sub></td><td>17.7</td><td>23.3</td></tr><tr><td>n</td><td>AQUAL global</td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>1.09<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0018<sup>+0.0009</sup><sub>-0.0010</sub></td><td>28.4<sup>+0.4</sup><sub>-0.4</sub></td><td>8.4</td><td>0.86<sup>+0.07</sup><sub>-0.06</sub></td><td>26.3</td><td>33.2</td></tr><tr><td>n</td><td>AQUAL local</td><td>1.19<sup>+0.06</sup><sub>-0.04</sub></td><td>1.31<sup>+0.05</sup><sub>-0.05</sub></td><td>0.0048<sup>+0.0094</sup><sub>-0.0020</sub></td><td>29.4<sup>+0.5</sup><sub>-0.5</sub></td><td>8.7</td><td>0.79<sup>+0.07</sup><sub>-0.06</sub></td><td>1100</td><td>1490</td></tr><tr><td>n</td><td>QUMOND global</td><td>1.03<sup>+0.04</sup><sub>-0.04</sub></td><td>1.09<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0049<sup>+0.0014</sup><sub>-0.0017</sub></td><td>28.4<sup>+0.4</sup><sub>-0.4</sub></td><td>8.4</td><td>0.86<sup>+0.07</sup><sub>-0.06</sub></td><td>20.9</td><td>25.6</td></tr><tr><td>n</td><td>QUMOND local</td><td>1.12<sup>+0.05</sup><sub>-0.05</sub></td><td>1.23<sup>+0.04</sup><sub>-0.04</sub></td><td>0.0051<sup>+0.0032</sup><sub>-0.0018</sub></td><td>29.1<sup>+0.4</sup><sub>-0.5</sub></td><td>8.6</td><td>0.82<sup>+0.07</sup><sub>-0.06</sub></td><td>1080</td><td>1470</td></tr></tbody></table>

dynamics should be detected if these models are correct. Note, however, that the expected $\alpha_{grav}$ is less than unity – the result for the Simple IF at $a_0 = 1.2 \times 10^{-10}$ m s$^{-2}$ – in all cases, largely due to the lower preferred $a_0$ value. The tensions with the quadrupole measurement of equation (2) are strong, ~8–9$\sigma$. The left panels of Figs 2 and 3 show the full PPDs of $Q_2$ and $\alpha_{grav}$, respectively.

### 4.1.2 Extended mass-to-light ratio priors

We begin our investigation of more general *M/L* models with a Gaussian hyperprior model in which the centres of the lognormal priors on $\Upsilon_{disc}$ and $\Upsilon_{bulge}$ are free parameters to be inferred alongside $a_0$ and shape. These results comprise the first two rows of Table 2 (for the $\delta$-family). We find the best-fitting values of the hyperprior centres to be $\mu_b \approx 0.7$ and 0.6 regardless of the IF family and EFE model used. This provides a significantly better fit than the fiducial SPARC model, and the BIC shows that the addition of these two extra parameters is clearly warranted by the data. These values are, however, highly unexpected from a stellar population modelling point of view: in dwarfs, one should if anything expect $\Upsilon_{disc} < 0.5$ (F. Lelli, private communication), while the older stars in the bulge should in all cases have higher *M/L* than those in the disc. In Fig. 4, we show partial corner plots of the Gaussian hyperprior model for the $\delta$-family and two different EFE models: no EFE (left panel) and AQUAL EFE with maximum-clustering prior (right panel).

In addition to improving the goodness of fit, allowing this extra freedom in $\Upsilon_{disc}$ and $\Upsilon_{bulge}$ significantly increases $\delta$ from ~1 to up to ~1.3–1.5, corresponding to a sharper deep-MOND-to-Newtonian transition. As a result, the predicted $Q_2$ and $\alpha_{grav}$ values are significantly reduced, as shown in the central panels of Figs 2 and 3. The effect is most pronounced using the AQUAL local EFE model. However, although the tensions with Cassini and the WB

of Banik et al. (2024) are eased relative to the fiducial *M/L* model they remain significant, at the ~7$\sigma$ level for $Q_2$.

The other *M/L* models of Section 3.2 – shown in the remaining rows of Table 2 – tell similar stories. Allowing $\Upsilon_{disc}$ and $\Upsilon_{bulge}$ to float freely yields $\Upsilon_{disc} > \Upsilon_{bulge}$ with a corresponding significant improvement in goodness-of-fit. The fact that BIC is significantly lower for the free uniform than free hyperprior model, along with the best-fitting scatter on $\Upsilon_{disc}$ and $\Upsilon_{bulge}$ between galaxies being considerably larger than 25 per cent (italicized entries), suggests that the *M/L* values may not in fact possess the degree of similarity expected in the SPARC error model, at least under the assumption of a fixed RAR. Forcing $\mu_d < \mu_b$ or $\Upsilon_{disc} < \Upsilon_{bulge}$ results in $\mu_d \approx \mu_b$ and $\Upsilon_{disc} \approx \Upsilon_{bulge}$, respectively, with worsened goodness-of-fit and increased $Q_2$ and $\alpha_{grav}$. We omit the $e_N$ constraints in this table because they are very similar to those in Table 1, and show only $\Delta$BIC [not $\Delta$BIC(P)] for goodness-of-fit.

Even more drastic changes are apparent in the final four rows of Table 2, where we remove either galaxies with bulges or galaxies without bulges (in both cases using a Gaussian hyperprior for the remaining *M/Ls*). Removing the galaxies with bulges causes an immense increase in best-fitting shape – reaching $\delta = 2.5$ for the AQUAL local EFE model with maximum-clustering prior – and a more modest decrease in $a_0$. The RAR points transformed according to the best-fitting galaxy parameter values with the best-fitting $\delta$-family fit overplotted is shown in Fig. 5, illustrating the very sharp transition between the Newtonian and deep-MOND regimes. The SS, with acceleration $2.32 \times 10^{-10}$ m s$^{-2}$ shown by the horizontal blue line, is almost fully Newtonian. (In contrast, the Simple IF with the same $a_0$, shown in red, gives a large MOND boost at that acceleration.) This reduces the expected $Q_2$ almost to zero as shown in the right panel of Fig. 2, bringing it approximately into consistency with the Cassini measurement. The expected $\alpha_{grav}$ is concomitantly reduced to 0 (right panel of Fig. 3), making it fully consistent with
<figure><img src="image_4.png" /><figcaption>**Figure 2.** PPD of $Q_2$ from the RAR for the $\delta$-family under the five different EFE models. Left panel: Fiducial SPARC M/L model; centre panel: Gaussian hyperprior model; right panel: as centre but excluding galaxies with bulges. The ‘local’ EFE models, with a separate $e_N$ per galaxy, use the maximum-clustering prior. The constraint on $Q_2$ from Cassini is visible in dashed grey in the right panel.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** As Fig. 2, but for $\alpha_{grav}$ rather than $Q_2$. The constraint on $\alpha_{grav}$ from the WBT of Banik et al. (2024) is visible in dashed grey in the right panel.</figcaption></figure>

the null result of Banik et al. (2024). The overall goodness-of-fit is, however, reduced, as shown by larger $\Delta$BIC and $\Delta$BIC(P) values. We note that the use of a Gaussian hyperprior for the mean of the disc M/L $\mu_d$ is important for the great increase in the best-fitting shape parameter when removing galaxies with bulges; there is almost no such increase if the fiducial SPARC M/L model is used instead, with best-fitting $\delta$ returning to $\sim 1$. Conversely, removing galaxies without bulges lowers shape, increases $a_0$, brings back a strong tension with the $Q_2$ measurement, and predicts $\alpha_{grav} \approx 1$.

These results indicate an incompatibility between the separate RAR fits of bulgey and bulge-free galaxies. Fig. 6 illustrates this by showing various fits for one bulgey galaxy, UGC 2953. The fits obtained when fixing the values of $a_0$ and shape parameter to their best fit values obtained from bulge-free galaxies ($\delta = 1.97$, $a_0 = 0.99 \times 10^{-10} \text{ m s}^{-2}$ without EFE, and $\delta = 2.49$, $a_0 = 1.1 \times 10^{-10} \text{ m s}^{-2}$ with EFE) are quite poor. In addition, the best-fitting distance, inclination and $\Upsilon$ values are far from their

priors. It therefore appears unlikely that the high shape values preferred by the bulge-free galaxies can explain the dynamics of bulgey galaxies. To provide another angle on this discrepancy, we show in Fig. 7 the {$a_0$, shape} constraints without EFE from the fiducial SPARC model, the free hyperprior model and the models with bulgey or bulge-free galaxies removed (using a Gaussian hyperprior for the remaining M/Ls). Each model is in clear tension with the others. This could be due either to the modelling of the bulge and disc light within the SPARC data-reduction pipeline, or to inaccuracy of the RAR modelling of these components. For example, the violation of axisymmetry caused by a significant bulge or bar – or the dependence on past trajectories and exact orbital structure in modified inertia formulations – may greatly modify the effective MOND force law. Using the fiducial M/L model rather than the Gaussian hyperprior for the no-bulge and bulge-only cases the best-fitting $\delta$ values are reduced, to $\sim 1$ for the former and $\sim 0.9$ for the latter, but they remain clearly dis-
<table><caption>**Table 2.** As Table 1, but focusing on the δ-family and models either without an EFE or with a galaxy-specific AQUAL EFE with maximum-clustering prior on $e_N$, allowing the M/L model to vary. The models considered are (1) unconstrained Gaussian hyperpriors on $\Upsilon_{disc}$ and $\Upsilon_{bulge}$, (2) as (1) but requiring the centre of the $\Upsilon_{bulge}$ prior to exceed that of the $\Upsilon_{disc}$ prior, (3) free uniform priors on $\Upsilon_{disc}$ and $\Upsilon_{bulge}$, (4) as (3) but requiring $\Upsilon_{bulge} > \Upsilon_{disc}$ galaxy-by-galaxy, (5) as (2) but excluding galaxies with bulges, and (6) as (2) but excluding galaxies without bulges. For models with Gaussian hyperpriors on $\Upsilon_{disc}$ and $\Upsilon_{bulge}$, the $\mu_d$ and $\mu_b$ columns show the mean and 68 per cent confidence intervals of their centres; for the other models (italicized) they instead show the mean and 68 per cent confidence intervals of the mean $\Upsilon_{disc}$ and $\Upsilon_{bulge}$ values across the posteriors of all galaxies. The uncertainty on $\sigma_{int}$ (measured in dex) is $\pm 0.001$ in all cases. $\Delta$BIC is again defined relative to the top row of Table 1.</caption><thead><tr><th>M/L model</th><th>EFE model</th><th>shape (δ)</th><th>a<sub>0</sub></th><th>σ<sub>int</sub></th><th>μ<sub>d</sub></th><th>μ<sub>b</sub></th><th>Q<sub>2</sub></th><th>σ<sub>Q2</sub></th><th>α<sub>grav</sub></th><th>ΔBIC</th></tr></thead><tbody><tr><td>Free hyper</td><td>No EFE</td><td>1.28<sup>+0.06</sup><sub>-0.06</sub></td><td>1.04<sup>+0.03</sup><sub>-0.03</sub></td><td>0.034</td><td>0.71<sup>+0.03</sup><sub>-0.03</sub></td><td>0.62<sup>+0.04</sup><sub>-0.03</sub></td><td>25.5<sup>+0.8</sup><sub>-0.9</sub></td><td>7.2</td><td>0.40<sup>+0.05</sup><sub>-0.05</sub></td><td>-53.0</td></tr><tr><td>Free hyper</td><td>AQUAL local</td><td>1.50<sup>+0.08</sup><sub>-0.07</sub></td><td>1.21<sup>+0.04</sup><sub>-0.04</sub></td><td>0.032</td><td>0.72<sup>+0.03</sup><sub>-0.03</sub></td><td>0.62<sup>+0.04</sup><sub>-0.03</sub></td><td>25.6<sup>+0.9</sup><sub>-0.9</sub></td><td>7.2</td><td>0.34<sup>+0.05</sup><sub>-0.05</sub></td><td>1000</td></tr><tr><td>Const hyper</td><td>No EFE</td><td>1.24<sup>+0.06</sup><sub>-0.06</sub></td><td>1.05<sup>+0.03</sup><sub>-0.03</sub></td><td>0.034</td><td>0.68<sup>+0.02</sup><sub>-0.02</sub></td><td>0.69<sup>+0.02</sup><sub>-0.02</sub></td><td>26.3<sup>+0.7</sup><sub>-0.8</sub></td><td>7.5</td><td>0.44<sup>+0.05</sup><sub>-0.05</sub></td><td>-46.9</td></tr><tr><td>Const hyper</td><td>AQUAL local</td><td>1.45<sup>+0.07</sup><sub>-0.07</sub></td><td>1.23<sup>+0.04</sup><sub>-0.04</sub></td><td>0.032</td><td>0.69<sup>+0.02</sup><sub>-0.02</sub></td><td>0.70<sup>+0.02</sup><sub>-0.02</sub></td><td>26.6<sup>+0.8</sup><sub>-0.9</sub></td><td>7.6</td><td>0.39<sup>+0.05</sup><sub>-0.05</sub></td><td>1000</td></tr><tr><td>Free unif</td><td>No EFE</td><td>1.28<sup>+0.06</sup><sub>-0.06</sub></td><td>1.10<sup>+0.03</sup><sub>-0.03</sub></td><td>0.031</td><td>0.73<sup>+0.72</sup><sub>-0.40</sub></td><td>0.68<sup>+0.41</sup><sub>-0.27</sub></td><td>26.5<sup>+0.7</sup><sub>-0.8</sub></td><td>7.6</td><td>0.44<sup>+0.05</sup><sub>-0.05</sub></td><td>-440</td></tr><tr><td>Free unif</td><td>AQUAL local</td><td>1.52<sup>+0.08</sup><sub>-0.07</sub></td><td>1.24<sup>+0.04</sup><sub>-0.04</sub></td><td>0.029</td><td>0.74<sup>+0.70</sup><sub>-0.38</sub></td><td>0.60<sup>+0.53</sup><sub>-0.19</sub></td><td>25.9<sup>+0.9</sup><sub>-0.9</sub></td><td>7.3</td><td>0.34<sup>+0.05</sup><sub>-0.05</sub></td><td>637</td></tr><tr><td>Const unif</td><td>No EFE</td><td>1.06<sup>+0.04</sup><sub>-0.04</sub></td><td>1.18<sup>+0.04</sup><sub>-0.04</sub></td><td>0.032</td><td>0.58<sup>+0.73</sup><sub>-0.33</sub></td><td>0.73<sup>+0.58</sup><sub>-0.27</sub></td><td>30.1<sup>+0.4</sup><sub>-0.5</sub></td><td>8.9</td><td>0.75<sup>+0.06</sup><sub>-0.06</sub></td><td>-274</td></tr><tr><td>Const unif</td><td>AQUAL local</td><td>1.24<sup>+0.06</sup><sub>-0.05</sub></td><td>1.39<sup>+0.05</sup><sub>-0.05</sub></td><td>0.031</td><td>0.57<sup>+0.72</sup><sub>-0.29</sub></td><td>0.74<sup>+0.63</sup><sub>-0.28</sub></td><td>30.8<sup>+0.5</sup><sub>-0.6</sub></td><td>9.1</td><td>0.68<sup>+0.06</sup><sub>-0.06</sub></td><td>823</td></tr><tr><td>No bulge</td><td>No EFE</td><td>1.97<sup>+0.22</sup><sub>-0.18</sub></td><td>0.99<sup>+0.03</sup><sub>-0.03</sub></td><td>0.041</td><td>0.81<sup>+0.03</sup><sub>-0.03</sub></td><td>—</td><td>14.3<sup>+3.0</sup><sub>-2.9</sub></td><td>2.7</td><td>0.05<sup>+0.04</sup><sub>-0.04</sub></td><td>2140</td></tr><tr><td>No bulge</td><td>AQUAL local</td><td>2.49<sup>+0.33</sup><sub>-0.27</sub></td><td>1.10<sup>+0.04</sup><sub>-0.04</sub></td><td>0.039</td><td>0.81<sup>+0.03</sup><sub>-0.03</sub></td><td>—</td><td>10.9<sup>+3.2</sup><sub>-2.9</sub></td><td>1.9</td><td>0.00<sup>+0.03</sup><sub>-0.02</sub></td><td>2940</td></tr><tr><td>Only bulge</td><td>No EFE</td><td>0.99<sup>+0.06</sup><sub>-0.06</sub></td><td>1.30<sup>+0.08</sup><sub>-0.07</sub></td><td>0.026</td><td>0.59<sup>+0.05</sup><sub>-0.05</sub></td><td>0.68<sup>+0.04</sup><sub>-0.04</sub></td><td>31.7<sup>+0.5</sup><sub>-0.6</sub></td><td>9.4</td><td>0.99<sup>+0.13</sup><sub>-0.12</sub></td><td>-191</td></tr><tr><td>Only bulge</td><td>AQUAL local</td><td>1.19<sup>+0.08</sup><sub>-0.07</sub></td><td>1.57<sup>+0.10</sup><sub>-0.09</sub></td><td>0.024</td><td>0.62<sup>+0.05</sup><sub>-0.04</sub></td><td>0.69<sup>+0.04</sup><sub>-0.04</sub></td><td>32.7<sup>+0.7</sup><sub>-0.8</sub></td><td>9.6</td><td>0.87<sup>+0.11</sup><sub>-0.10</sub></td><td>-105</td></tr></tbody></table>

<table><caption>**Table 3.** As Table 2, but using the modified gravity prescription of equation (18) for calculating the MONDian radial acceleration $g$ rather than the algebraic MOND relation. Here, we use the *n*-family without EFE in all cases.</caption><thead><tr><th>M/L model</th><th>shape (n)</th><th>a<sub>0</sub></th><th>σ<sub>int</sub></th><th>μ<sub>d</sub></th><th>μ<sub>b</sub></th><th>Q<sub>2</sub></th><th>σ<sub>Q2</sub></th><th>α<sub>grav</sub></th><th>ΔBIC</th></tr></thead><tbody><tr><td>Fiducial</td><td>0.78<sup>+0.03</sup><sub>-0.03</sub></td><td>0.96<sup>+0.05</sup><sub>-0.05</sub></td><td>0.046</td><td>0.5</td><td>0.7</td><td>29.2<sup>+0.3</sup><sub>-0.2</sub></td><td>8.7</td><td>1.34<sup>+0.09</sup><sub>-0.08</sub></td><td>2050</td></tr><tr><td>Free hyper</td><td>1.57<sup>+0.09</sup><sub>-0.09</sub></td><td>1.12<sup>+0.03</sup><sub>-0.03</sub></td><td>0.042</td><td>0.97<sup>+0.03</sup><sub>-0.03</sub></td><td>0.64<sup>+0.04</sup><sub>-0.04</sub></td><td>22.8<sup>+1.0</sup><sub>-1.0</sub></td><td>6.3</td><td>0.33<sup>+0.05</sup><sub>-0.05</sub></td><td>2290</td></tr><tr><td>Free unif</td><td>1.34<sup>+0.08</sup><sub>-0.07</sub></td><td>1.10<sup>+0.04</sup><sub>-0.03</sub></td><td>0.041</td><td>0.88<sup>+0.76</sup><sub>-0.40</sub></td><td>0.64<sup>+0.43</sup><sub>-0.25</sub></td><td>25.2<sup>+0.8</sup><sub>-0.9</sub></td><td>7.1</td><td>0.49<sup>+0.06</sup><sub>-0.06</sub></td><td>2120</td></tr><tr><td>No bulge</td><td>1.93<sup>+0.19</sup><sub>-0.17</sub></td><td>1.01<sup>+0.03</sup><sub>-0.03</sub></td><td>0.049</td><td>1.08<sup>+0.04</sup><sub>-0.04</sub></td><td>—</td><td>17.5<sup>+2.4</sup><sub>-2.4</sub></td><td>3.8</td><td>0.15<sup>+0.05</sup><sub>-0.04</sub></td><td>3680</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Partial corner plots of the RAR inference using the δ-family of IFs and the Gaussian hyperprior model for $\Upsilon_{disc}$ and $\Upsilon_{bulge}$, for the case of no EFE (left panel) and AQUAL with maximum-clustering $e_N$ prior (right panel). The truth line shows $\mu_b = 0.7$, the fiducial SPARC value (the corresponding $\mu_d = 0.5$ is off the plot).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** SPARC data points transformed according to the best-fitting parameters for the δ-family, local AQUAL EFE, and M/L hyperprior model removing galaxies with bulges. The orange line and band show the best-fitting RAR and its 95 per cent confidence limit (producing a very small boost to gravity in the SS), while the red line shows the Simple IF with the same $a_0$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Example of fits for an individual bulgey galaxy, UGC 2953. The blue curves correspond to a fit where $a_0$ and shape are fixed to their best-fitting values based on the RAR of bulgey galaxies only (last lines of Table 2). The orange curves correspond instead to a fit where $a_0$ and shape are fixed to their best-fitting values based on the RAR of bulgeless galaxies only. The solid lines correspond to fits including the AQUAL-local EFE, while the dashed lines correspond to fits without EFE. The fitted parameters (distances, inclination, and $\gamma$) are additionally pushed far from their priors for the orange fits. This illustrates that bulgey galaxies are not well described by the RAR obtained from bulge-free galaxies alone, which could otherwise satisfy the $Q_2$ and WBT constraints.</figcaption></figure>

crepant in $a_0$. We leave further investigation of this issue to further work.

### 4.1.3 AQUAL modified gravity RAR fit

Our fiducial analysis above uses the algebraic MOND relation, equation (3), which is also the usual form of the RAR. While this is an exact prediction in modified inertia formulations of the theory for circular orbits, modified gravity formulations necessarily deviate from it to a (small) extent that depends on the geometry of the system (see Famaey & McGaugh 2012). As we are interested here in testing modified gravity theories (specifically AQUAL and QUMOND) in which the IF is the same between galaxy RCs, the SS, and nearby wide binaries, it is important to check that our results are not significantly

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of constraints on $a_0$ and $\delta$ from the RAR under various M/L models, assuming no EFE. In the context of MOND, the discordance between the ellipses of the fiducial, Gaussian hyperprior, no-bulge, and bulge-only fits seem to hint at systematics in the SPARC data, subtleties related to the underlying MOND theory and/or significant issues with the fiducial M/L priors. This disagreement cannot be resolved by changing the EFE model.</figcaption></figure>

affected by this deviation from the algebraic relation. To this end, we employ here the AQUAL approximation found by Brada & Milgrom (1995, equation 25) for flat discs:

$$g = \frac{g_N}{\mu \left( \frac{g_N^+}{a_0} \nu \left( \frac{g_N^+}{a_0} \right) \right)}, \quad (18)$$

where $g_N^+ = (g_N^2 + (2\pi G \Sigma)^2)^{1/2}$, for a total baryonic surface density $\Sigma$ at each point in the disc. Here, $\mu$ is the usual alternative form of the IF already described in Section 3.4, defined by $\nu(x\mu(x)) = 1/\mu(x)$. equation (18) is appropriate for the radial acceleration predicted by AQUAL within razor-thin axisymmetric systems; we caution that its accuracy will be lower for galaxies with bulges.

The baryonic surface density is the sum of a stellar disc, stellar bulge and gas component; while the first two are given in the SPARC catalogue, the third is not. We calculate it here from the SPARC $V_{gas}$ values using equation (13) of Toomre (1963). This requires an integral of $V_{gas}$ to $r = \infty$. We extrapolate it beyond the last measured point with a Keplerian decline, corresponding to the assumption of spherical symmetry and that all gas is enclosed within the observed region. To test the effect of this approximation, we consider an alternative (extreme) model in which $V_{gas}$ is constant beyond the last measured point, finding minimal difference in the results.

We use the MONDian $g$ calculated from equation (18) in place of equation (3) in the likelihood term and repeat the inference for the $n$-family (which has simple analytic forms for both $\nu$ and $\mu$), without EFE. The results are shown in Table 3. We see a decrease in the quality of the fit (shown by $\sigma_{int}$ and $\Delta$BIC exceeding the corresponding rows of Table 1 and 2) and an increased disagreement between the RAR, $Q_2$ and WBT constraints (shown by the reduced shape and increased $\sigma_{Q_2}$ and $\alpha_{grav}$). Our main results using the algebraic MOND model are therefore conservative. It is also interesting to note that the RAR seems to favour the straight algebraic MOND relation over this modified gravity correction, although not with high significance.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Constraints on $a_0$, $\delta$, and $g_{\text{ext}}$ from the Cassini constraint $Q_2 = (3 \pm 3) \times 10^{-27} \text{ m s}^{-2}$, using a flat prior on $Q_2$. $a_0$ and $g_{\text{ext}}$ are in units of $10^{-10} \text{ m s}^{-2}$. The $g_{\text{ext}}$ posterior is dominated by its truncated Gaussian prior.</figcaption></figure>

## 4.2 Constraints from Cassini and comparison to the RAR and WBT

To shed further light on the discrepancy at the heart of our analysis, we now infer {$a_0$, shape} from the SS quadrupole. To set the stage, we show in Fig. 8 the corner plot of the $Q_2$ inference for the $\delta$-family when $g_{\text{ext}}$ is also inferred using a truncated Gaussian prior from Gaia Collaboration et al. (2021) (see Section 3.3). Higher shape and lower $a_0$ values are preferred by the likelihood because this pushes the SS further into the Newtonian regime where the prediction for $Q_2$ is 0. However, our flat prior on $Q_2$ creates a prior $\propto d(Q_2)/d(\text{shape})$ on shape, which goes to 0 at high shape as $Q_2$ levels out at 0. This truncates the posterior so that only a finite shape range is allowed. The shape and $a_0$ marginals only decline at small values (also in Fig. 9) due to our priors $\delta > 0.5$ and $a_0 > 0.0025 \times 10^{-10} \text{ m s}^{-2}$, which do not impact the (in)consistency with the RAR. There is little constraining power on $g_{\text{ext}}$ from the $Q_2$ measurement, so its posterior is similar to its prior.

Next, in Fig. 9, we compare the posteriors on $a_0$ and shape from the RAR and $Q_2$ inferences. In each panel, we use the $\delta$-family and show the results for each of the EFE models; the three panels show the same *M/L* models as in Figs 2 and 3. We see again the clear tension between the RAR and $Q_2$ measurement using either the fiducial or Gaussian hyperprior *M/L* model: The inferences only become consistent when the galaxies with bulges are removed. This conclusion is not affected by the choice of IF family or EFE model. If instead of sampling $g_{\text{ext}}$ we choose the value in the range $2-2.48 \times 10^{-10} \text{ m s}^{-2}$ that minimizes $Q_2$ at given $a_0$ and shape (corresponding to the best possible $g_{\text{ext}}$ for reducing the tension), the posterior is only slightly shifted to lower shape.

It is also instructive to consider the RAR fits obtained when the $Q_2$ constraint is enforced at the $1\sigma$ level. To achieve this we fix $\delta = 5$ and $a_0 = 1.45 \times 10^{-10} \text{ m s}^{-2}$ – which yields $Q_2 = 6 \times 10^{-27} \text{ s}^{-2}$ – and repeat the RAR inference using Gaussian hyperpriors on $\Upsilon_{\text{disc}}$ and $\Upsilon_{\text{bulge}}$. For the no-EFE model, requiring this IF lowers $\ln(\hat{\mathcal{L}})$ by 195, indicating a hugely worse fit. The $\ln$-distance prior is also worsened by 16.4 and the $\ln$-inclination prior by 38.8, showing that

these galaxy parameters are being forced to take values far from their priors. The constraints on the *M/L* hyperpriors are $\mu_d = 0.92 \pm 0.03$ and $\mu_b = 0.67 \pm 0.04$, indicating an even greater deviation from the fiducial *M/L* model than the regular Gaussian hyperprior case. $\sigma_{\text{int}}$ is only 0.04 dex, however, and the plot of the transformed points with $\delta = 5$, $a_0 = 1.45 \times 10^{-10} \text{ m s}^{-2}$ model overlaid (not shown) looks similar to Fig. 5, just with a somewhat sharper transition. Analogous results hold for the other EFE models.

Finally, we show in Fig. 10 the PPD of $\alpha_{\text{grav}}$ from the $Q_2$ chain. Here, we see that the Cassini measurement already imposes a stronger constraint on $\alpha_{\text{grav}}$ than does the WBT result of Banik et al. (2024): Equation (2) allows one to predict, from the SS alone, that the WBT will be null to high precision regardless of IF family. This result may appear surprising given the fact that Banik et al. rule out the RAR IF with $a_0 = 1.2 \times 10^{-10} \text{ m s}^{-2}$ at $16\sigma$, while we rule out the RAR IF using the Cassini measurement at only $8.7\sigma$ (top row of Table 1). This is primarily because $\alpha_{\text{grav}}$ is more sensitive to shape than is $Q_2$, so when allowing shape to increase from unity one achieves $\alpha_{\text{grav}} \approx 0$ sooner than $Q_2 \approx 0$. We therefore conclude that as a general constraint on modified gravity MOND, the SS quadrupole is stronger than the WBT that is currently possible. The PPD of $\alpha_{\text{grav}}$ from $Q_2$ agrees with the inference of Banik et al. but would not with those of Hernandez et al. (2019, 2022, 2023), Hernandez (2023), and Chae (2023a, 2023b), who find approximate agreement with the Simple or RAR IF, $\alpha_{\text{grav}} \approx 1$. Our results therefore argue in favour of Banik et al.

## 4.3 Can a fine-tuned IF reconcile the Cassini and RAR constraints?

One might be tempted to design a new IF with a sharp transition between the MONDian and Newtonian regimes in an attempt to reconcile the Cassini and RAR constraints, whilst still producing a deviation from Newtonian gravity in local wide binaries. We have seen that removing bulges and adopting free hyperpriors for the stellar *M/Ls* of SPARC discs can reconcile the RAR and the Cassini constraints to better than $2\sigma$. However, this produces no non-Newtonian behaviour in the WBT. Moreover, using fiducial SPARC *M/L* values or including bulge galaxies in the RAR fit returns us to a gradual transition close to the Simple or RAR IF.

In order to explore this question more systematically, let us now consider a class of IFs built from $v_{\text{RAR}}$ by introducing a new acceleration scale $a_{\text{N, trans}}$ where a sharp sigmoid transition is applied between $v_{\text{RAR}}$ and the Newtonian regime ($v = 1$). This means that for $y < a_{\text{N, trans}}/a_0$ the IF behaves as $v_{\text{RAR}}$, but it becomes fully Newtonian as soon as $y > a_{\text{N, trans}}/a_0$. The case $a_{\text{N, trans}} = a_0$ would correspond to a sharp transition from Newton to deep-MOND without introducing a new acceleration scale. However, to keep a positive WBT, one would need $a_{\text{N, trans}} > a_0$ in order to produce non-Newtonian behaviour at the Sun's position in the Milky Way.

The left panel of Fig. 11 shows the RAR obtained with this class of fine-tuned IFs for different values of $a_{\text{N, trans}}$. We also display the minimum and maximum observed accelerations probed by the SPARC RCs (using the cuts defined in Section 2.1), and the external field of the Milky Way at the position of the sun (assuming $a_0 = 1.03 \times 10^{-10} \text{ m s}^{-2}$; see Table 1). The right panel then shows the value of $Q_2$ in the SS as a function of $a_{\text{N, trans}}$. It is clear that, in order to satisfy the Cassini constraint at $1\sigma$, the new transition acceleration $a_{\text{N, trans}}$ must be below the Newtonian external field at the Sun. This will both strongly affect the RAR fits, especially in bulge galaxies that do not allow such a sharp transition as we have seen in Section 4.1.2, and necessarily yield no deviation from Newton in
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Posteriors on $a_0$ and $\delta$ from the RAR and SS quadrupole. Left panel: fiducial SPARC M/L model; centre panel: free Gaussian hyperpriors on $\Upsilon_{\text{disc}}$ and $\Upsilon_{\text{bulge}}$; right panel: as centre but excluding galaxies with bulges.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** PPD of $\alpha_{\text{grav}}$ for the three IF families from the $Q_2$ constraint. This constraint foretells $\alpha_{\text{grav}} = 0$ to higher precision than the WBT of Banik et al. (2024).</figcaption></figure>

local wide binaries. It appears therefore to be highly unlikely that one could cook up a fine-tuned universal IF that would simultaneously explain the RAR and pass the Cassini constraint, especially if one also wants a non-Newtonian WBT.

## 5 DISCUSSION AND CONCLUSION

The simplicity and regularity of galaxy kinematics – epitomized for late-type galaxies’ radial dynamics by the RAR – appears to provide strong support for MOND relative to the galaxy formation scenario of ΛCDM. We find excellent fits to the RAR of the SPARC galaxy sample using generalized IFs of the n-, δ-, or γ-families, with the very small intrinsic scatter (~8 per cent) approximately expected from remaining systematics or departures from spherical symmetry in modified gravity formulations of MOND. Including all galaxies, we find the best-fitting shape of the RAR to be very close to the Simple IF or the one proposed in McGaugh et al. (2016). The impact of the EFE, while not clearly required in our analysis, is consistent with the expectation from independent calculations of the galaxies’ large-scale environments.

It is interesting to compare our best-fitting $a_0$ values with those in the literature. We recover $a_0 \approx 1 \times 10^{-10} \text{ m s}^{-2}$ when not including

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Exploration of a class of IFs constructed by introducing a sharp sigmoid transition between $v_{\text{RAR}}$ and $v = 1$ (fully Newtonian behaviour) at a new transition acceleration scale $a_{\text{N, trans}}$. Left panel: the RAR of this IF class for different values of $a_{\text{N, trans}}$. The two horizontal dotted lines bound the gravitational fields probed by the SPARC RAR while the dashed horizontal line shows the external field acting on the SS, $g_{\text{ext, Sun}}$. The vertical shaded area shows the range of Newtonian external fields corresponding to $g_{\text{ext, Sun}}$ for different IFs. Right panel: the variation of $Q_2$ produced by this class of IFs as a function of $a_{\text{N, trans}}$. This transition must occur at an acceleration below the external field $g_{\text{N, ext, Sun}}$ in order to satisfy the Cassini constraint at $1\sigma$, which would both strongly impact the SPARC RAR and produce no non-Newtonian behaviour in the WBT (since local wide binaries are embedded in the same external field as the SS, where the IF would have to be fully Newtonian).</figcaption></figure>

the EFE or when the EFE is fitted as a global parameter, and $a_0 \approx 1.2 \times 10^{-10} \text{ m s}^{-2}$ when the EFE is allowed to vary galaxy-by-galaxy with the maximum clustering prior. Historically, the first estimates of $a_0$ – fitted with the Standard IF until circa 2005 – gave values around $a_0 = 1.2 \times 10^{-10} \text{ m s}^{-2}$ (e.g. Kent 1987; Milgrom 1988; Begeman, Broeils & Sanders 1991), typically using fewer than 10 RCs. Complementing the data of Begeman et al. (1991) with data from Gentile et al. (2004) and using the Simple IF, Famaey et al. (2007a) found $a_0 = 1.35 \times 10^{-10} \text{ m s}^{-2}$ with 14 RCs, whilst a later analysis of 12 higher resolution RCs from the THINGS survey (Gentile et al. 2011) gave $a_0 = 1.22 \times 10^{-10} \text{ m s}^{-2}$ with the Simple IF, but with an intriguing tension for the archetypal RC of NGC 3198, indicating $a_0 = 0.9 \times 10^{-10} \text{ m s}^{-2}$ in line with the more recent EFE-free estimates. Including the local EFE gives back the ‘old’ value of $a_0 \approx 1.2 \times 10^{-10} \text{ m s}^{-2}$.

Our main result is that the location and sharpness of the transition between the deep-MOND and Newtonian regimes inferred from the
RAR is grossly inconsistent with that inferred from the Cassini measurement of the SS quadrupole in classical modified gravity formulations of MOND. Including full marginalization over all relevant galaxy nuisance parameters, we calculate this tension to be 8.7σ for the δ-family of IFs and fiducial M/L model of SPARC. We use the algebraic MOND relation between $g_{\text{obs}}$ and $g_{\text{bar}}$ fiducially, but show explicitly that switching to an AQUAL fitting formula does not qualitatively affect the results. It is also interesting to note that the RAR seems to slightly favour the straight algebraic MOND relation over this modified gravity correction, although not with high significance.

We find that the tension can be slightly ameliorated by allowing the M/L $Υ$ of the discs and bulges of the SPARC galaxies to float, which is strongly favoured by the RAR data (given an IF-family fit). This significantly increases the sharpness of the MOND transition – making the Solar neighbourhood more Newtonian – but prefers $Υ_{\text{disc}} > Υ_{\text{bulge}}$, flying in the face of stellar population modelling. Using Gaussian hyperpriors with free means for $Υ_{\text{disc}}$ and $Υ_{\text{bulge}}$ reduces the $Q_2$ tension to $7.2σ$. A much greater gain is to be had by discarding the 31 galaxies with bulges, in which case the RAR with free $Υ_{\text{disc}}$ is in only $1.9σ$ tension with the quadrupole constraint using the AQUAL galaxy-by-galaxy EFE model. The resulting model cannot fit galaxies with bulges, however. We stress that these modifications to the quadrupole prediction arise entirely internally to the RAR modelling and without any reference to the Cassini measurement: It is not a case of ameliorating tension in a joint inference by enlarging the parameter space, but rather a preference within the SPARC data. We also show that, for the wide range of parametric RAR shapes considered here, the Cassini measurement implies (in the modified gravity context) no detectable deviation from Newtonian gravity in wide binaries in the Solar neighbourhood to greater precision than Gaia measurements of the wide binaries currently afford. This is consistent with (and renders unsurprising) the results of Banik et al. (2024), but disagrees with Hernandez et al. (2019, 2022, 2023), Hernandez (2023), Chae (2023a, 2023b), and Hernandez & Chae (2023).

Thus, while the RAR *does* appear to be tantamount to a natural law in galaxies (Lelli et al. 2017; Desmond 2023; Stiskalek & Desmond 2023), extending it into the Solar neighbourhood does not work unless the SPARC M/L model and/or treatment of bulges is greatly in error. How should we interpret this perplexing situation? Aside from the conclusion that the simplicity of galaxy dynamics is a red herring – a to-be-understood emergent behaviour of dark matter – and MOND on a hiding to nothing, there appear to be three possibilities:

(i) Our different results from considering separately the bulge and bulge-free galaxies – and preference for $Υ$s far removed from the fiducial SPARC values – may be indicating that there are systematics in the RAR data which could permit a MOND-to-Newtonian transition as steep as is required by the SS quadrupole. Removing galaxies with bulges and allowing $Υ_{\text{disc}}$ to float freely is the most conservative model if one has reason to doubt the bulge modelling, and it may not be a coincidence that it is consistent with the Cassini measurement. In this case, the orbits of WBs in the Solar neighbourhood are also expected to be fully Newtonian. This may, however, cause problems with observables not studied here, such as the RC of the Milky Way.

There is no a priori reason to remove galaxies with bulges, but a posteriori the tension between galaxies with and without bulges provides strong motivation to examine this issue further (within the MOND paradigm). The fact that the RAR $\{a_0, \text{shape}\}$ constraints between bulgey and bulge-free galaxies and for various M/L choices

are clearly discrepant with one another (Fig. 7) is concerning in itself for the universality of the RAR. One may wonder whether the RAR parameters may correlate with other galaxy properties; this appears not to be the case from Stiskalek & Desmond (2023). Note that the bulgey versus bulge-free dichotomy is different to that of Chae (2022), where it is instead between the RARs followed by inner and outer points of the RCs regardless of bulge fraction. This is unlikely to remain in the underlying RARs presented here due to their extremely small intrinsic scatter.

(ii) Remaining within the modified gravity paradigm and taking the inconsistency of the RAR and $Q_2$ results at face value necessitates a modification to the AQUAL or QUMOND weak-field limits of any such theory. One way to achieve this would be to screen the modification to GR on small scales, a method already common for preventing scalar-tensor theories of gravity from manifesting unobserved deviations from the inverse square law ('fifth forces') in the SS (see Jain & Khoury 2010; Baker et al. 2021 for reviews). This may be achieved in the case of MOND by, for instance, including a covariant Galileon-type term (involving a length-scale) for the scalar field in the MOND action (Babichev, Deffayet & Esposito-Farèse 2011): This could completely suppress all MOND effects below a Vainshtein radius depending on the source mass, $a_0$ and the new length-scale. Interestingly, with the length-scale proposed by Babichev et al. (2011), MOND effects would be killed in the SS and up to the typical (mass and length) scale of WBs. Such a Galileon term could therefore be added to frameworks such as those of Skordis & Złośnik (2021), and would not affect the predicted RAR on larger scales where interesting deviations from the algebraic MOND relation are already to be found. This is studied for high-mass (galaxy cluster-like) systems in Durakovic & Skordis (2023); see also Banik et al. (2024). The idea of adding a new scale to MOND is somewhat akin to the approach of 'Extended MOND' where a second gravitational variable (in this case potential) modulates the MOND boost (Zhao & Famaey 2012).

It is important to bear in mind that our conclusions refer solely to specific modified gravity instantiations of MOND in the classical regime, namely the AQUAL and QUMOND models where the IF is necessarily universal. Generalizations such as TRIMOND (Milgrom 2023b) and GQUMOND (Milgrom 2023c) have more complex behaviour, and may evade our constraints by having different IFs in galaxies, wide binaries and the SS. In GQUMOND, a new scale naturally arises, allowing MOND effects to be screened below some length-scale (as in the Galileon screening case) or dynamical time. TRIMOND manifests three gravitational degrees of freedom – the MOND potential and two auxiliary potentials, one of which is the Newtonian potential – and contains AQUAL and QUMOND as special cases (Milgrom 2023b). There thus remains an extensive theory space of MOND-like theories, which could be constrained through the apparent incompatibility revealed here. However, one may consider that we are reaching the point where the complexity of modified gravity formulations of MOND is no longer warranted by its empirical successes relative to ΛCDM.

(iii) Finally, one may consider MOND as implying a modification to inertia. Modified inertia formulations are far more difficult to construct and study because the dynamics of an object may depend on its entire past history (e.g. Milgrom 2011; Milgrom 2022). It is, however, modified inertia that produces in its pristine form the strongest piece of MOND phenomenology, viz. the RAR (Milgrom 1994; Milgrom 2022). In this context, the mismatch between bulgey and bulge-free galaxies could be a manifestation of the subtle dependence of the MOND force law on the underlying orbital dynamics and history. Our results may therefore be indicating that modified inertia is the
correct interpretation of MOND, spurring the search for theories and models that allow it to make contact with data beyond disc galaxy dynamics. Indeed, an SS quadrupole is not generically predicted by modified inertia Milgrom (2023a), potentially allowing the Cassini constraint to be completely circumvented.

## ACKNOWLEDGEMENTS

We thank Indranil Banik, Michal Bilek, Amel Durakovic, Federico Lelli, Stacy McGaugh, Mordehai Milgrom, Constantinos Skordis, Richard Stiskalek and Tariq Yasin for useful discussions.

HD is supported by a Royal Society University Research Fellowship (grant no. 211046). BF acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 834148) and from the Agence Nationale de la Recherche (ANR projects ANR-18-CE31-0006 and ANR-19-CE31-0017).

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 693024).

For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

## DATA AVAILABILITY

All data underlying this paper may be shared on request to the authors.

## REFERENCES

* Adelberger E. G., Heckel B. R., Nelson A. E., 2003, Annu. Rev. Nucl. Part. Sci., 53, 77

* Angus G. W., Famaey B., Buote D. A., 2008, MNRAS, 387, 1470

* Antreasian P. G. et al., 2005, AAS/AIAA Astrodynamics Conference. Jet Propulsion Laboratory, National Aeronautics and Space Administration, Pasadena, CA, p. AAS 05-312

* Asencio E., Banik I., Mieske S., Venhola A., Kroupa P., Zhao H., 2022, MNRAS, 515, 2981

* Babichev E., Deffayet C., Esposito-Farèse G., 2011, Phys. Rev. D, 84, 061502

* Baker T. et al., 2021, Rev. Mod. Phys., 93, 015003

* Banik I., Zhao H., 2022, Symmetry, 14, 1331

* Banik I., Thies I., Famaey B., Candlish G., Kroupa P., Ibata R., 2020, ApJ, 905, 135

* Banik I., Pittordis C., Sutherland W., Famaey B., Ibata R., Mieske S., Zhao H., 2024, MNRAS, 527, 4573

* Begeman K. G., Broeils A. H., Sanders R. H., 1991, MNRAS, 249, 523

* Bekenstein J., Milgrom M., 1984, ApJ, 286, 7

* Bingham E. et al., 2019, J. Mach. Learn. Res., 20, 28:1

* Blanchet L., Novak J., 2011, MNRAS, 412, 2530

* Bovy J., Bird J. C., García Pérez A. E., Majewski S. R., Nidever D. L., Zasowski G., 2015, ApJ, 800, 83

* Brada R., Milgrom M., 1995, MNRAS, 276, 453

* Brown K., Mathur H., 2023, AJ, 166, 168

* Chae K.-H., 2022, ApJ, 941, 55

* Chae K.-H., 2023a, ApJ, 960, 114

* Chae K.-H., 2023b, ApJ, 952, 128

* Chae K.-H., Milgrom M., 2022, ApJ, 928, 24

* Chae K.-H., Lelli F., Desmond H., McGaugh S. S., Li P., Schombert J. M., 2020, ApJ, 904, 51

* Chae K.-H., Desmond H., Lelli F., McGaugh S. S., Schombert J. M., 2021, ApJ, 921, 104

* Chae K.-H., Lelli F., Desmond H., McGaugh S. S., Schombert J. M., 2022, Phys. Rev. D, 106, 103025

* de Blok W. J. G., McGaugh S. S., 1997, MNRAS, 290, 533

* Desmond H., 2017a, MNRAS, 464, 4160

* Desmond H., 2017b, MNRAS, 472, L35

* Desmond H., 2023, MNRAS, 526, 3342

* Desmond H., Wechsler R. H., 2015, MNRAS, 454, 322

* Desmond H., Katz H., Lelli F., McGaugh S., 2019, MNRAS, 484, 239

* Desmond H., Bartlett D. J., Ferreira P. G., 2023, MNRAS, 521, 1817

* Di Cintio A., Lelli F., 2016, MNRAS, 456, L127

* Durakovic A., Skordis C., 2023, preprint (arXiv:2312.00889)

* Famaey B., Binney J., 2005, MNRAS, 363, 603

* Famaey B., McGaugh S. S., 2012, Living Rev. Relativ., 15, 10

* Famaey B., Gentile G., Bruneton J.-P., Zhao H., 2007a, Phys. Rev. D, 75, 063002

* Famaey B., Bruneton J.-P., Zhao H., 2007b, MNRAS, 377, L79

* Famaey B., McGaugh S., Milgrom M., 2018, MNRAS, 480, 473

* Folkner W. M., Williams J. G., Boggs D. H., Park R., Kuchynka P., 2014, IPN Progress Report 42

* Freundlich J., Famaey B., Oria P.-A., Bilek M., Müller O., Ibata R., 2022, A&A, 658, A26

* Gaia Collaboration et al., 2021, A&A, 649, A9

* Gentile G., Salucci P., Klein U., Vergani D., Kalberla P., 2004, MNRAS, 351, 903

* Gentile G., Famaey B., de Blok W. J. G., 2011, A&A, 527, A76

* Ghari A., Haghi H., Zonoozi A. H., 2019a, MNRAS, 487, 2148

* Ghari A., Famaey B., Laporte C., Haghi H., 2019b, A&A, 623, A123

* Haghi H., Baumgardt H., Kroupa P., Grebel E. K., Hilker M., Jordi K., 2009, MNRAS, 395, 1549

* Haghi H., Bazkiaei A. E., Zonoozi A. H., Kroupa P., 2016, MNRAS, 458, 4172

* Haghi H. et al., 2019, MNRAS, 487, 2441

* Hees A. et al., 2012, Class. Quant. Grav., 29, 235027

* Hees A., Folkner W. M., Jacobson R. A., Park R. S., 2014, Phys. Rev. D, 89, 102002

* Hees A., Famaey B., Angus G. W., Gentile G., 2016, MNRAS, 455, 449

* Hernandez X., 2023, MNRAS, 525, 1401

* Hernandez X., Chae K.-H., 2023, preprint (arXiv:2312.03162)

* Hernandez X., Cortés R. A. M., Allen C., Scarpa R., 2019, Int. J. Mod. Phys. D, 28, 1950101

* Hernandez X., Cookson S., Cortés R. A. M., 2022, MNRAS, 509, 2304

* Hernandez X., Verteletskyi V., Nasser L., Aguayo-Ortiz A., 2023, MNRAS, 528, 4720

* Iess L., Rappaport N. J., Jacobson R. A., Racioppa P., Stevenson D. J., Tortora P., Armstrong J. W., Asmar S. W., 2010, Science, 327, 1367

* Iess L. et al., 2012, Science, 337, 457

* Jacobson R. A. et al., 2006, AJ, 132, 2520

* Jain B., Khoury J., 2010, Ann. Phys., 325, 1479

* Jeffreys H., 1939, The Theory of Probability, Oxford Classic Texts in the Physical Sciences, Oxford, England

* Jones-Smith K., Abraham R., Kell L., Mathur H., 2018, New Journal of Physics, 20, 063042, https://dx.doi.org/10.1088/1367-2630/aaca23

* Kent S. M., 1987, AJ, 93, 816

* Kroupa P. et al., 2018, Nature, 561, E4

* Kroupa P. et al., 2022, MNRAS, 517, 3613

* Kuzio de Naray R., Kaufmann T., 2011, MNRAS, 414, 3617

* Lelli F., McGaugh S. S., Schombert J. M., 2016a, AJ, 152, 157

* Lelli F., McGaugh S. S., Schombert J. M., 2016b, ApJ, 816, L14

* Lelli F., McGaugh S. S., Schombert J. M., Pawlowski M. S., 2017, ApJ, 836, 152

* Li P. et al., 2023, A&A, 677, A24

* McGaugh S., Milgrom M., 2013, ApJ, 775, 139

* McGaugh S. S., Schombert J. M., Bothun G. D., de Blok W. J. G., 2000, ApJ, 533, L99

* McGaugh S. S., Lelli F., Schombert J. M., 2016, Phys. Rev. Lett., 117, 201101

* McGaugh S. S., Lelli F., Schombert J. M., 2020, Res. Notes AAS, 4, 45

* Maquet L., Pierret F., 2015, Phys. Rev. D, 91, 084015

* Milgrom M., 1983a, ApJ, 270, 365

* Milgrom M., 1983b, ApJ, 270, 371

* Milgrom M., 1983c, ApJ, 270, 384

* Milgrom M., 1988, ApJ, 333, 689

* Milgrom M., 1994, Ann. Phys., 229, 384
Milgrom M., 2009, MNRAS, 399, 474

Milgrom M., 2010, MNRAS, 403, 886

Milgrom M., 2011, Acta Phys. Polon. B, 42, 2175

Milgrom M., 2022, Phys. Rev. D, 106, 064060

Milgrom M., 2023a, preprint (arXiv:2310.14334)

Milgrom M., 2023b, Phys. Rev. D, 108, 063009

Milgrom M., 2023c, Phys. Rev. D, 108, 084005

Oman K. A. et al., 2015, MNRAS, 452, 3650

Oria P. A. et al., 2021, ApJ, 923, 68

Paranjape A., Sheth R. K., 2021, MNRAS, 507, 632

Phan D., Pradhan N., Jankowiak M., 2019, preprint (arXiv:1912.11554)

Pittordis C., Sutherland W., 2023, The Open J. Astrophys., 6, 4

Pointecouteau E., Silk J., 2005, MNRAS, 364, 654

Roper F. A., Oman K. A., Frenk C. S., Benítez-Llambay A., Navarro J. F., Santos-Santos I. M. E., 2023, MNRAS, 521, 1316

Sanders R. H., 1999, ApJ, 512, L23

Schombert J., McGaugh S., 2014, Publ. Astron. Soc. Aust., 31, e036

Schwarz G., 1978, Ann. Stat., 6, 461

Skordis C., Złośnik T., 2021, Phys. Rev. Lett., 127, 161302

Stiskalek R., Desmond H., 2023, MNRAS, 525, 6130

Swaters R. A., Sancisi R., van Albada T. S., van der Hulst J. M., 2009, A&A, 493, 871

The L. S., White S. D. M., 1988, AJ, 95, 1642

Thomas G. F., Famaey B., Ibata R., Renaud F., Martin N. F., Kroupa P., 2018, A&A, 609, A44

Toomre A., 1963, ApJ, 138, 385

Tully R. B., Fisher J. R., 1977, A&A, 54, 661

Vokrouhlický D., Nesvorný D., Tremaine S., 2024, preprint (arXiv:2403.09555)

Wolf A., Smith J., 1995, Control Eng. Pract., 3, 1611

Wu X., Famaey B., Gentile G., Perets H., Zhao H., 2008, MNRAS, 386, 2199

Zhao H. S., Famaey B., 2006, ApJ, 638, L9

Zhao H., Famaey B., 2012, Phys. Rev. D, 86, 067301

Zwaan M. A., van der Hulst J. M., de Blok W. J. G., McGaugh S. S., 1995, MNRAS, 273, L35

This paper has been typeset from a *TEX/LATEX* file prepared by the author.