On the tension between the Radial Acceleration Relation and Solar System quadrupole in modified gravity MOND

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 AQUAL and 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 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.


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 et al. 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 ≲ 10 −10 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 ⋆ harry.desmond@port.ac.uk † aurelien.hees@obspm.fr‡ benoit.famaey@astro.unistra.frexample, 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 et al. 2016b;Desmond et al. 2019), which one would typically expect to be produced by the dark matter halos of a Λ-Cold Dark Matter (Λ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 Λ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 Λ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 noncircular 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 et al. 2016;Lelli et al. 2017) between the total gravitational acceleration at any point in the disc (g obs ) and that generated by the observed baryons (g 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,b,c) promoted them to a tentative law of Nature by proposing that Newtonian dynamics, the weak-field limit of GR, breaks down at ultra-low accelerations typical of the outskirts of galaxies (below a new universal constant a0 ≃ 10 −10 m/s 2 ) where one instead has g obs = √ a0 g bar .In semi-analytic and hydrodynamical models of galaxy formation of Λ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 et al. 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 a0, and of bodies in the Solar System (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 a0 (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 a0.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(Hernandez et al. , 2022;;Hernandez 2023;Hernandez et al. 2023;Chae 2024Chae , 2023;;Hernandez & Chae 2023).
In the present work, we focus instead on the SS tests, which were revolutionised by the Cassini spacecraft's radioscience tracking (Wolf & Smith 1995;Antreasian et al. 2005;Jacobson et al. 2006;Iess et al. 2010Iess et al. , 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 MONDthe 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 et al. 2007b;Wu et al. 2008;Haghi et al. 2009;McGaugh & Milgrom 2013;Hees et al. 2016;Haghi et al. 2016;Famaey et al. 2018;Kroupa et al. 2018;Thomas et al. 2018;Haghi et al. 2019;Chae et al. 2020a;Banik et al. 2020;Oria et al. 2021;Kroupa et al. 2022;Asencio 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;Milgrom 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 obs = f (g 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 bar ) → g bar as g bar → ∞ and f (g bar ) → √ a0 g bar as g bar → 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 et al. ( 2023) evaluated all possible functions of low complexity on dynamical galaxy data, concluding that the optimal function for RC data is fairly complexand 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 obs ∝ √ g bar at low g bar .Ad-ditional uncertainty arises from the EFE, which depends on the imperfectly observationally characterised baryonic largescale 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 a0 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 massto-light ratios 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.Second, we compare the RAR {a0, 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 Sec. 2 we describe the galaxy and SS data that we use as constraints.In Sec. 3 we lay out our MOND models (Sec.3.1), our method for inferring their parameters from the RAR (Sec.3.2) and SS quadrupole (Sec.3.3), and the method for relating these to the WBT (Sec.3.4).Sec. 4 describes the results, separately for the RAR (Sec.4.1) and SS quadrupole (Sec.4.2), and explores why the two sets of constraints are difficult to reconcile (Sec.4.3).We conclude in Sec. 5. Throughout the paper, log has base 10.

SPARC galaxy sample
For the RAR we utilise the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample (Lelli et al. 2016a), 1 containing 175 RCs from the literature with photometry at 3.6µm from the Spitzer satellite (see also Gentile et al. 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 Hi distributions) or inclination i < 30 • , 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).

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) where ê = gext/gext is a unit vector pointing towards the Galactic Centre, x the position within the SS with respect to the Sun, δij the Kronecker delta and Q2 a parameter that depends on the MOND IF and acceleration scale a0, the Sun's mass M and the value of the external gravitational field from the Galaxy, gext.
The modified Newtonian potential from Eq. 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ý et al. 2024).In Hees et al. (2014), the Q2 parameter of Eq. 1 was inferred using the DE430 planetary ephemerides data (Folkner et al. 2014).Of central importance for constraining this parameter are the 9 years 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σ constraint may be described by (2)

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: where g is the total gravitational field (or dynamical acceleration), g N is the Newtonian gravitational field sourced by baryons, ν(y) is the IF and gN ≡ |g N |.The only a priori requirement on ν from the basic tenets of MOND is the asymptotic behaviour ν(y) → 1 when y → ∞, and ν(y) → 1/ √ y when y → 0. The three most common examples are the "Simple" IF (Famaey & Binney 2005;Zhao & Famaey 2006) the "Standard" IF (Milgrom 1983c): (5) and the "RAR" (also called "McGaugh-Lelli-Schombert" or "MLS") IF (McGaugh et al. 2016) 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 gN/a0, their only other degree of freedom) that interpolates between specific IFs by controlling the sharpness of the transition between the deep-MOND (gN/a0 ≪ 1) and Newtonian (gN/a0 ≫ 1) regimes.Here, we investigate three IF families described in Famaey & McGaugh (2012): 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 δ-and γ-families contain the RAR IF at δ = γ = 1.2 We refer to n, δ and γ 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" (Freundlich et al. 2022;Oria et al. 2021) and "AQUAL numerical" models of Chae & Milgrom (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 a0, i.e.Milgrom 2022).These roughly span the space of possible EFE behaviour, although note that Eq. 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 eN of the SPARC galaxies.The first is simply to treat eN as a global constant and infer it (with a wide uniform prior) from the data.The second is to allow the local eN to vary galaxy-by-galaxy, with a prior specified by the "maximum clustering" model of Chae et al. (2021), which infers eN for each SPARC galaxy from its large-scale baryonic environment assuming unseen baryons to correlate maximally with observed baryons.3eN for each galaxy separately is then a free parameter of the inference to be marginalised 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 eN 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(eN) = −2.300± 0.575.
Finally, it is important to note that in MOND Eq. 3 (often referred to as the "algebraic MOND relation") is exact only for circular orbits in modified inertia versions of the paradigm (Milgrom 1994;Milgrom 2022), and cannot be strictly exact in modified gravity, even in the absence of an EFE.We will consider it a sufficient approximation (Brown et al. 2018;Oria et al. 2021) in our fiducial analysis, but in Sec.4.1.3we will investigate an alternative formula for the MOND boost taking into account the correction for flattened systems within the AQUAL formulation.

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 (a0 and intrinsic scatter σ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 L3.6, the massto-light ratios (M/L) of the disc Υ disk , bulge Υ bulge and gas Υgas, and the external field strengths eN.This amounts to ∼900 parameters in total, which is too high a dimensionality for most Markov Chain Monte Carlo (MCMC) 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 (NUTS, a species of Hamiltonian Monte Carlo) as implemented in numpyro (Phan et al. 2019;Bingham et al. 2019).This produces a fully converged chain of ∼2000 points in ∼10 minutes; we concatenate 28 chains with multiprocessing for improved statistics and to verify insensitivity to the initialisation.Restricting to the Simple IF and a single parametrisation of the EFE, Desmond (2023) revealed the intrinsic scatter of the RAR to be minute (σint = 0.034 ± 0.001(stat) ± 0.001(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 mass-to-light ratios of the SPARC galaxies' disks 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 con-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. (2021Chae et al. ( , 2022)), and is a priori the most likely case in MOND.The differences for our purposes are minor.
straining the RAR parameters.In particular we consider six possible priors: (i) The fiducial SPARC model in which Υ disk and Υ bulge follow lognormal4 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) Υ disk and Υ bulge drawn from separate Gaussian hyperpriors with means µ d and µ 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 µ b > µ d , as expected from population synthesis models.(iv) Υ disk and Υ bulge drawn independently from wide uniform priors, modelling a scenario in which nothing is known about them a priori.(v) As (iv) but requiring Υ bulge > Υ disk 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 bulgey and bulge-free galaxies.5 The Υ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 et al. 2019a).
We use uniform priors on a0, σ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 a0 or σint, or a prior on shape corresponding to a flat prior on Q2 (see Sec. 3.3) makes negligible difference to the RAR results.

Quadrupole inference
Theoretically, the SS quadrupole value that appears in Eq. 1 can be expressed as (Milgrom 2009) where, in the SS, M = 1 M⊙, and q is a dimensionless parameter that depends only on the MOND IF as well as the value of where gext 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: where ν = ν e 2 N + v 4 + 2eNv 2 ξ and eN is the solution of eNν (eN) = ẽ.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 andBlanchet &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 Q2,6 and is therefore conservative when it comes to assessing the tension with the data.
The value of q(ẽ) obtained from Eq. 12 depends mainly on the behaviour of the IF around eN .To illustrate this, Fig. 1 presents, as a function of the argument of the ν-function in a small bin between y0 and y0 + ∆y (with ∆y set to 0.1), the contribution ∆q(y0) of ν(y) in this bin to the value of q(ẽ).In practice, ∆q(y0) is computed by replacing the function ν(y) by 1 in Eq. 12 (i.e.g = gN) except in the range [y0, y0 + ∆y] where it instead follows Eq. 6 (i.e., δ = 1 in Eq. 7b or γ = 1 in Eq. 7c).Fig. 1 presents the variation of ∆q with y0 for three different values of the external field ẽ.The vertical dashed lines show the corresponding values of eN .It can be noticed that |∆q| is large only around eN , which is the case for all other IFs too.This means that the Cassini constraint probes the IF around y ∼ eN = gN,ext/a0, 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 gN,ext in these different cases.
Then, for each IF family (Eq.7), we have computed the q factor on a regular grid for the (ẽ, shape) parameters with ẽ 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 a0, we supplement this with an irregular sampling of ẽ values extending to 1000.We then interpolate this grid using a 2D cubic spline.We use Eq. 2 as the likelihood function in our inference to constrain {a0, 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 gext = 2.32 ± 0.16 × 10 −10 m s −2 (Gaia Collaboration et al. 2021).While the 2σ lower bound 2 × 10 −10 m s −2 is realistic (corresponding to a large but acceptable peculiar velocity of the Sun; Bovy et al. 2015), the 2σ upper bound is not (as it would imply an unrealistic negative peculiar velocity of the Sun).We therefore consider the upper bound on gext to be 2.48 × 10 −10 m s −2 .We treat it in one of two ways: i) use the value within the range [2, 2.48]×10 −10 m s −2 that leads to the lowest predicted value of Q2 in order to be conservative in our inference of {a0, shape} and maximise consistency with the RAR, or ii) infer gext along with {a0, shape} using the prior N (2.32, 0.16) truncated at 2.48.These methods give very similar results; we present results for the

e=2
Figure 1.Variation of ∆q(y 0 ), the contribution from ν(y) between y 0 and y 0 +0.1 to q(ẽ), for the IF ν RAR of Eq. 6 and three different values of ẽ.The value of q(ẽ) from Eq. 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(ẽ) mainly probes the behaviour of the IF around e N .
second method because it is more statistically principled and makes the Q2 and RAR inferences fully independent.
For the Q2 inference, a flat prior on shape produces a posterior extending to infinity, in which limit Q2 → 0. As well as complicating the construction of confidence intervals, this would introduce a strong volume effect into the posterior predictive distribution (PPD) of the WBT test statistic (see next subsection).We therefore instead adopt a prior flat on Q2 (with a uniform prior on a0), which matches that used in Hees et al. (2014) to derive the constraint of Eq. 2. This is achieved by taking the numerical derivative of Q2 with respect to shape at each point within the grid.

Connection to the wide binary test
Following Banik et al. (2024), we summarise the expected result of a WBT in the Solar neighbourhood using a parameter αgrav defined as where η ≡ ⟨gr⟩/gN, the boost of the azimuthally averaged asymptotic radial gravity compared to Newton.In QUMOND this is given by where νe ≡ ν (eN) is the IF evaluated at the Newtonianequivalent Galactic gravity eN = gN,ext/a0 at the Solar position.In AQUAL it is instead given by µ is an alternative description of the IF, defined by µ(x) = ν(y) −1 where y = xµ(x).Eq. 13 is chosen such that αgrav = 0 for fully Newtonian gravity, and αgrav = 1 for the case of QUMOND with the Simple IF, assuming a0 = 1.2 × 10 −10 m s −2 and ẽ = 1.8.Table 4 of Banik et al. (2024) shows the αgrav values corresponding to a few choices of IF at fixed a0 and gN,ext.
Eq. 13 describes a mapping from {a0, shape, gN,ext} to αgrav, allowing us to convert results on these parameters from either the RAR or Q2 into a PPD of αgrav values.This is the distribution we would expect to see given the RAR or Q2 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 ge = 2.142 × 10 −10 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 {a0, shape} posteriors from the inference of Banik et al., as we do in Sec.4.2 from the Q2 measurement.This is in order not to restrict ourselves to those results, reflecting the less established nature of the WBT than the Q2 measurement and the fact that other WBT analyses have reached substantially different conclusions (Hernandez et al. 2019(Hernandez et al. , 2022;;Hernandez 2023;Hernandez et al. 2023;Chae 2023Chae , 2024)).
Like Q2, αgrav goes to 0 in the limit shape → ∞, so a flat prior on shape and the associated semi-infinite prior volume towards that limit would have made αgrav appear more strongly constrained by Cassini than it in fact is.This is cured by adopting a flat prior on Q2.Note that the αgrav inference of Banik et al. (2024), to which we compare the αgrav PPD from Q2, used instead a uniform prior on αgrav.We have checked that switching to this prior-or a log-uniform prior in Q2 or αgrav-makes little difference to our results.

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 σ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 eN varies per galaxy, we show the median and 68% confidence intervals of the median eN values across the posteriors of all galaxies; these entries are italicised to distinguish them from the qualitatively different global eN constraints appearing in the same column.For each model we also summarise the PPDs of Q2 and αgrav by their median and 68 per cent confidence interval, including the tension with the Cassini measurement in the case of Q2.This is calculated assuming a Gaussian distribution of the lower uncertainty.
The final two columns describe the goodness-of-fit, specifically the 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  and 9).For the models with galaxy-specific ("local") EFE the quoted e N constraints-written in italics-describe the median stacked posteriors over all galaxies, implicitly including the maximum-clustering prior.The Q 2 (Eq.10) and αgrav (Eq.13) columns summarise the posterior predictive distributions of the quadrupole and WB test statistic from the RAR fits, while σ Q2 is the number of sigma tension between the Q 2 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.
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) where k is the number of free parameters, N the number of data points and L the maximum likelihood value.
In the limit of much data, the likelihood approximates the posterior because the prior does not scale with the number of datapoints 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 for maximum posterior value 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: • a0 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 δ-and γ-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 δ-family because it is the simplest that includes the most popular RAR function.Some historical context to our a0 constraints may be found in Sec. 5.  1) but requiring the centre of the Υ bulge prior to exceed that of the Υ disk prior, 3) free uniform priors on Υ disk and Υ bulge , 4) as (3) but requiring Υ bulge > Υ disk 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 Υ disk and Υ bulge , the µ d and µ b columns show the mean and 68% confidence intervals of their centres; for the other models (italicised) they instead show the mean and 68% confidence intervals of the mean Υ disk and Υ bulge values across the posteriors of all galaxies.The uncertainty on σ int (measured in dex) is ±0.001 in all cases.
∆BIC is again defined relative to the top row of Table 1. that while the AQUAL constraints on a global eN are significantly below the expectation from the maximum-clustering prior (and hence also the galaxy-by-galaxy eN 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. (2020bChae et al. ( , 2021Chae et al. ( , 2022) ) due to their focus on specific galaxies and ours on the overall RAR, and our different fitting and goodness-of-fit assessment procedures.
• Q2 and αgrav are clearly non-zero in all cases, indicating that deviation from Newtonian gravity in the SS quadrupole and WB dynamics should be detected if these models are correct.Note however that the expected αgrav is less than unity-the result for the Simple IF at a0 = 1.2 × 10 −10 m s −2 -in all cases, largely due to the lower preferred a0 value.
The tensions with the quadrupole measurement of Eq. 2 are strong, ∼8-9σ.The left panels of Figs. 2 and 3 show the full PPDs of Q2 and αgrav respectively.

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 Υ disk and Υ bulge are free parameters to be inferred alongside a0 and shape.These results comprise the first two rows of Table 2 (for the δ-family).We find the best-fit values of the hyperprior centres to be µ d ≈ 0.7 and µ b ≈ 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 Υ disk < 0.5 (F.Lelli, priv. comm.), while the older stars in the bulge should in all cases have higher M/L than those in the disk.In Fig. 4 we show partial corner plots of the Gaussian hyperprior model for the δ-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 Υ disk and Υ bulge significantly increases δ from ∼1 to up to ∼1.3−1.5, corresponding to a sharper deep-MOND-to-Newtonian transition.As a result, the predicted Q2 and α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 WBT of Banik et al. (2024) are eased relative to the fiducial M/L model they remain significant, at the ∼7σ level for Q2.
The other M/L models of Sec.3.2-shown in the remaining rows of Table 2-tell similar stories.Allowing Υ disk and Υ bulge to float freely yields Υ disk > Υ 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-fit scatter on Υ disk and Υ bulge between galaxies being considerably larger than 25 per cent (italicised 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 µ d < µ b or Υ disk < Υ bulge results in µ d ≈ µ b and Υ disk ≈ Υ bulge respectively, with worsened goodness-offit and increased Q2 and αgrav.We omit the eN constraints in this table because they are very similar to those in Table 1, and show only ∆BIC (not ∆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-fit shapereaching δ = 2.5 for the AQUAL local EFE model with maximum-clustering prior-and a more modest decrease in a0.The RAR points transformed according to the best-fit galaxy parameter values with the best-fit δ-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 × 10 −10 m s −2 shown by the horizontal blue line, is almost fully Newtonian.(In contrast, the Simple IF with the same a0, shown in red, gives a large MOND boost at that acceleration.)This reduces the expected Q2 almost to zero as shown in the right panel of Fig. 2, bringing it approximately into consistency with the Cassini measurement.The expected αgrav is concomitantly reduced to 0 (right panel of Fig. 3), making it fully consistent with the null result of Banik et al. (2024).The overall goodness-of-fit is however reduced, as shown by larger ∆BIC and ∆BIC(P) values.We note that the use of a Gaussian hyperprior for the mean of the disc mass-to-light ratio µ d is important for the great increase in the best-fit 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-fit δ returning to ∼1.Conversely, removing galaxies without bulges lowers shape, increases a0, brings back a strong tension with the Q2 measurement, and predicts αgrav ≈ 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 a0 and shape parameter to their best fit values obtained from bulge-free galaxies (δ = 1.97, a0 = 0.99 × 10 −10 m s −2 without EFE, and δ = 2.49, a0 = 1.1 × 10 −10 m s −2 with EFE) are quite poor.In addition, the best-fit distance, inclination and Υ values are far from their priors.It therefore appears unlikely that the high shape values preferred by the bulgefree galaxies can explain the dynamics of bulgey galaxies.To provide another angle on this discrepancy, we show in Fig. 7 the {a0, 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-fit δ values are reduced, to ∼1 for the former and ∼0.9 for the latter, but they remain clearly discrepant in a0.We leave further investigation of this issue to further work.

AQUAL modified gravity RAR fit
Our fiducial analysis above uses the algebraic MOND relation, Eq. 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 Solar System, and nearby wide binaries, it is important to check that our results are not significantly affected by this deviation from the algebraic relation.To this end, we employ here the AQUAL approximation found by Brada & Milgrom (1995, eq 25) for flat disks: where g + N ≡ g 2 N + (2πGΣ) 2 1/2 , for a total baryonic surface density Σ at each point in the disk.Here µ is the usual alternative form of the IF already described in Sect.3.4, defined by ν(xµ(x)) ≡ 1/µ(x).Eq. 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 disk, 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 Vgas values using eq 13 of Toomre (1963).This requires an integral of Vgas to r = ∞.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 Vgas is constant beyond the last measured point, finding minimal difference in the results.
We use the MONDian g calculated from Eq. 18 in place of Eq. 3 in the likelihood term and repeat the inference for the n-family (which has simple analytic forms for both ν and µ), without EFE.The results are shown in Table 3.We see a decrease in the quality of the fit (shown by σint and ∆BIC exceeding the corresponding rows of Table 1 and 2) and an increased disagreement between the RAR, Q2 and WBT constraints (shown by the reduced shape and increased σQ 2 and α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.

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 {a0, shape} from the SS quadrupole.To set the stage, we show in Fig. 8 the corner plot of the Q2 inference for the δ-family when gext is also inferred using a truncated Gaussian prior from Gaia Collaboration et al. ( 2021) (see Sec. 3.3).Higher shape and lower a0 values are preferred by the likelihood because this pushes the SS further into the  Newtonian regime where the prediction for Q2 is 0. However, our flat prior on Q2 creates a prior ∝ d(Q2)/d(shape) on shape, which goes to 0 at high shape as Q2 levels out at 0. This truncates the posterior so that only a finite shape range is allowed.The shape and a0 marginals only decline at small values (also in Fig. 9) due to our priors δ > 0.5 and a0 > 0.0025 × 10 −10 m s −2 , which do not impact the (in)consistency with the RAR.There is little constraining power on gext from the Q2 measurement, so its posterior is similar to its prior.Next, in Fig. 9 we compare the posteriors on a0 and shape from the RAR and Q2 inferences.In each panel we use the δ-family and show the results for each of the EFE models; the three panels show the same M/L models as in Figs. 2  2).The orange curves correspond instead to a fit where a 0 and shape are fixed to their best-fit 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 Υs) 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.and 3. We see again the clear tension between the RAR and Q2 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 Comparison of constraints on a 0 and δ 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.
instead of sampling gext we choose the value in the range 2 − 2.48 × 10 −10 m s −2 that minimises Q2 at given a0 and shape (corresponding to the best possible gext 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 Q2 constraint is enforced at the 1σ level.To achieve this we fix δ = 5 and a0 = 1.45 × 10 −10 m s −2 -which yields Q2 = 6 × 10 −27 s −2 -and repeat the RAR inference using Gaussian hyperpriors on Υ disk and Υ bulge .For the no-EFE model, requiring this IF lowers ln( 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 µ d = 0.92 ± 0.03 and µ b = 0.67 ± 0.04, indicating an even greater deviation from the fiducial M/L model than the regular Gaussian hyperprior case.σint is only 0.04 dex, however, and the plot of the transformed points with δ = 5, a0 = 1.45 × 10 −10 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 αgrav from the Q2 chain.Here we see that the Cassini measurement already imposes a stronger constraint on αgrav than does the WBT result of Banik et al. (2024): Eq. 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 a0 = 1.2 × 10 −10 m s −2 at 16σ, while we rule out the RAR IF using the Cassini measurement at only 8.7σ (top row of Table 1).This is primarily because αgrav is more sensitive to shape than is Q2, so when allowing shape to increase from unity one achieves αgrav ≈ 0 sooner than Q2 ≈ 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 αgrav from Q2 agrees with the inference of Banik et al. but  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/L ratios of SPARC disks can reconcile the RAR and the Cassini constraints to better than 2σ.However, this produces no non-Newtonian behaviour in the WBT.Moreover, using fiducial SPARC M/L values or including bulgey 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 νRAR by introducing a new acceleration scale aN,trans where a sharp sigmoid transition is applied between νRAR and the Newtonian regime (ν = 1).This means that for y < aN,trans/a0 the IF behaves as νRAR, but it becomes fully Newtonian as soon as y > aN,trans/a0.The case aN,trans = a0 would correspond to a sharp transition from Newton to deep-MOND without introducing a new acceleration scale.However, to keep a pos- Figure 10.PPD of αgrav for the three IF families from the Q 2 constraint.This constraint foretells αgrav = 0 to higher precision than the WBT of Banik et al. (2024).
itive WBT, one would need aN,trans > a0 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 aN,trans.We also display the minimum and maximum observed accelerations probed by the SPARC RCs (using the cuts defined in Sec.2.1), and the external field of the Milky Way at the position of the sun (assuming a0 = 1.03 × 10 −10 m s −2 ; see Table 1).The right panel then shows the value of Q2 in the SS as a function of aN,trans.It is clear that, in order to satisfy the Cassini constraint at 1σ, the new transition acceleration aN,trans must be below the Newtonian external field at the Sun.This will both strongly affect the RAR fits, especially in bulgey galaxies that do not allow such a sharp transition as we have seen in Sec.4.1.2,and necessarily yield no deviation from Newton in 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.

DISCUSSION AND CONCLUSION
The simplicity and regularity of galaxy kinematicsepitomised 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 generalised IFs of the n-, δ-or γ-families, with the very small intrinsic scatter (∼8 per cent) approximately expected from remaining systematics or departures from spherical symme-try in modified gravity formulations of MOND.Including all galaxies, we find the best-fit 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 external field effect, 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-fit a0 values with those in the literature.We recover a0 ≈ 1 × 10 −10 m s −2 when not including the EFE or when the EFE is fitted as a global parameter, and a0 ≈ 1.2 × 10 −10 m s −2 when the EFE is allowed to vary galaxy-by-galaxy with the maximum clustering prior.Historically, the first estimates of a0-fitted with the Standard IF until circa 2005-gave values around a0 = 1.2 × 10 −10 m s −2 (e.g.Kent 1987;Milgrom 1988;Begeman et al. 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 a0 = 1.35 × 10 −10 m s −2 with 14 RCs, whilst a later analysis of 12 higher resolution RCs from the THINGS survey (Gentile et al. 2011) gave a0 = 1.22×10 −10 m s −2 with the Simple IF, but with an intriguing tension for the archetypal RC of NGC 3198, indicating a0 = 0.9×10 −10 m s −2 in line with the more recent EFE-free estimates.Including the local EFE gives back the "old" value of a0 ≈ 1.2 × 10 −10 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 marginalisation 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 obs and g 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 mass-to-light ratios Υ of the disks 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 Υ disk > Υ bulge , flying in the face of stellar population modelling.Using Gaussian hyperpriors with free means for Υ disk and Υ bulge reduces the Q2 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 Υ disk 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. (2019Hernandez et al. ( , 2022)); Hernandez (2023); Hernandez et al. (2023);Chae (2024Chae ( , 2023)); 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 bulgey 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 Υ disk 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 {a0, 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 Q2 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 et al. 2011): this could completely suppress all MOND effects below a Vainshtein radius depending on the source mass, a0 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.Generalisations 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 2011Milgrom , 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, a Solar System quadrupole is not generically predicted by modified inertia Milgrom (2023a), potentially allowing the Cassini constraint to be completely circumvented.

Figure 2 .
Figure 2. PPD of Q 2 from the RAR for the δ-family under the five different EFE models.Left: Fiducial SPARC M/L model; centre: Gaussian hyperprior model; right: 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.

Figure 3 .
Figure 3.As Fig. 2, but for αgrav rather than Q 2 .The constraint on αgrav from the WBT of Banik et al. (2024) is visible in dashed grey in the right panel.

Figure 4 .
Figure 4. Partial corner plots of the RAR inference using the δ-family of IFs and the Gaussian hyperprior model for Υ disk and Υ bulge , for the case of no EFE (left) and AQUAL with maximum-clustering e N prior (right).The truth line shows µ b = 0.7, the fiducial SPARC value (the corresponding µ d = 0.5 is off the plot).

Figure 5 .
Figure5.SPARC data points transformed according to the bestfit parameters for the δ-family, local AQUAL EFE and M/L hyperprior model removing galaxies with bulges.The orange line and band show the best-fit 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 .
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-fit values based on the RAR of bulgey galaxies only (last lines of Table2).The orange curves correspond instead to a fit where a 0 and shape are fixed to their best-fit 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 Υs) 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.

Figure
Figure7.Comparison of constraints on a 0 and δ 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.
would not with those of Hernandez et al. (2019, 2022); Hernandez (2023); Hernandez et al. (2023); Chae (2024, 2023), who find approximate agreement with the Simple or RAR IF, αgrav ≈ 1.Our results therefore argue in favour of Banik et al.. 4.3 Can a fine-tuned IF reconcile the Cassini and RAR constraints?

Figure 9 .
Figure 9. Posteriors on a 0 and δ from the RAR and SS quadrupole.Left: Fiducial SPARC M/L model; centre: free Gaussian hyperpriors on Υ disk and Υ bulge ; right: as centre but excluding galaxies with bulges.

Figure 11 .
Figure 11.Exploration of a class of IFs constructed by introducing a sharp sigmoid transition between ν RAR and ν = 1 (fully Newtonian behaviour) at a new transition acceleration scale a N,trans .Left: the RAR of this IF class for different values of a 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 ext,Sun .The vertical shaded area shows the range of Newtonian external fields corresponding to g ext,Sun for different IFs.Right: the variation of Q 2 produced by this class of IFs as a function of a N,trans .This transition must occur at an acceleration below the external field g N,ext,Sun in order to satisfy the Cassini constraint at 1σ, 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).

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 0 has units of 10 −10 m s −2 and Q 2 of 10 −27 s −2 .We show all combinations of IF family (Eq.7) and EFE model (Eqs.8

Table 2 .
As Table1, 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 Υ disk and Υ bulge , 2) as ( • 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 eN values of the global-EFE models are all consistent with 0 within 2σ.It is however curious to note

Table 3 .
As Table2, but using the modified gravity prescription of Eq. 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.