The underlying radial acceleration relation

The radial acceleration relation (RAR) of late-type galaxies relates their dynamical acceleration, $g_\text{obs}$, to that sourced by baryons alone, $g_\text{bar}$, across their rotation curves. Literature fits to the RAR have fixed the galaxy parameters on which the relation depends -- distance, inclination, luminosity and mass-to-light ratios -- to their maximum a priori values with an uncorrelated Gaussian contribution to the uncertainties on $g_\text{bar}$ and $g_\text{obs}$. In reality these are free parameters of the fit, contributing systematic rather than statistical error. Assuming a range of possible functional forms for the relation with or without intrinsic scatter (motivated by Modified Newtonian Dynamics with or without the external field effect), I use Hamiltonian Monte Carlo to perform the full joint inference of RAR and galaxy parameters for the Spitzer Photometry and Accurate Rotation Curves (SPARC) dataset. This reveals the intrinsic RAR underlying that observed. I find an acceleration scale $a_0=(1.19 \pm 0.04 \, \text{(stat)} \pm 0.09 \, \text{(sys)}) \: \times \: 10^{-10}$ m s$^{-2}$, an intrinsic scatter $\sigma_\text{int}=(0.034 \pm 0.01 \, \text{(stat)} \pm 0.01 \, \text{(sys)})$ dex (assuming the SPARC error model is reliable) and weak evidence for the external field effect. I make summary statistics of all my analyses publicly available for future SPARC studies or applications of a calibrated RAR, for example redshift-independent distance measurement.


INTRODUCTION
Galaxies are observed to follow several tight and regular scaling relations between their internal motions and morphology.The classical correlations are the Tully-Fisher relation between rotation velocity and mass or luminosity in late-type galaxies (e.g.Tully & Fisher 1977;McGaugh et al. 2000;Pizagno et al. 2007) and the Fundamental Plane relating luminosity, size and velocity dispersion (e.g.Djorgovski & Davis 1987;Dressler et al. 1987;Cappellari et al. 2013, including its projection onto the mass-velocity plane, the Faber-Jackson relation; Faber & Jackson 1976) in early types.These are largely subsumed in late-type galaxies by the mass discrepancyacceleration or radial acceleration relation (Milgrom 1983a;Sanders 1990;McGaugh 2004;Lelli et al. 2017), relating the local total acceleration, g obs , to that sourced by baryons, g bar , across rotation curves.This provides more detailed radial information about the gravitational potential.
These relations provide the key evidence concerning the mass discrepancy problem in galaxies, namely that the motions of stars and gas imply far higher dynamical than baryonic masses in a Newtonian analysis.In the prevailing Λ Cold Dark Matter (ΛCDM) cosmology this difference is assumed to be made up by dark matter, leading to attempts to explain the relations through the modelling of galaxy forma-⋆ E-mail: harry.desmond@port.ac.uk tion, the galaxy-halo connection and halo mass distributions (e.g.Gnedin et al. 2007;Blanton et al. 2008;Desmond & Wechsler 2015, 2017;Di Cintio & Lelli 2016;Ludlow et al. 2017;Navarro et al. 2017;Keller & Wadsley 2017;Desmond 2017;Tenneti et al. 2018;Paranjape & Sheth 2021).However, the fact that galaxy formation in ΛCDM proceeds in a highly stochastic and complicated manner may make it difficult to explain "simple" (power-law or roughly double power-law) dynamical scaling relations.
An alternative hypothesis is that Newtonian gravity breaks down at the galaxy scale.Surprisingly, galaxy dynamics can be explained well by a model in which g obs = g bar for g bar ≫ a0 and g obs ∝ g 1/2 bar as g bar ≪ a0, where a0 ≈ 10 −10 m s −2 is a new fundamental constant.This naturally leads to the observed simplicity in the aforementioned scaling relations.Supplemented by an "interpolating function" that connects the Newtonian and modified gravity regimes, this theory is known as Modified Newtonian Dynamics (MOND ;Milgrom 1983a,c,b) and has achieved some success at explaining and even predicting galaxy behaviour (e.g.Famaey & McGaugh 2012;McGaugh & Milgrom 2013;Chae et al. 2020b).MOND has been incorporated into a range of nonrelativistic and relativistic theories over the past four decades, as reviewed most recently in Banik & Zhao (2022).
The RAR is MOND written in terms of observables for late-type galaxies.This at once gives the relation central importance in the missing mass debate and makes it the most sensitive probe of gravitational parameters within the MOND paradigm.These include the acceleration constant a0 marking the onset of modified dynamics, the intrinsic scatter σint and possibly a parameter eN ≡ gext/a0 describing the influence of mass surrounding the galaxy (external field effect, EFE ;Milgrom 1983a), where gext is the strength of the gravitational field in which the galaxy is embedded.σint bears on the question of whether the RAR manifests law-like gravitational behaviour as posited by MOND, and also determines the precision with which the relation may be used to calibrate galaxy properties such as distance (analogously to the Tully-Fisher relation).eN addresses the key question within the MOND paradigm of the extent to which-and manner in which-modified gravity or inertia violates the strong equivalence principle.Previous fits have found a0 ≈ 1.2 × 10 −10 m/s 2 , eN ≈ 0.003 and a small intrinsic scatter σint < 0.1 dex (Lelli et al. 2017;Li et al. 2018;Chae et al. 2021Chae et al. , 2022)).The existence of the EFE is however by no means well-established (for example Hernandez et al. 2019 andFreundlich et al. 2022 find evidence against it), and qualitatively similar phenomenology may arise in ΛCDM (Paranjape & Sheth 2022).
g bar and g obs depend on a number of properties of the galaxies, most importantly their distance D, inclination i, luminosity L and mass-to-light ratios Υ of their various components.These are nuisance parameters when determining the properties of the RAR, although of course of interest in their own right.Past RAR studies have either fixed these to their maximum a priori values given other measurements and then propagated their uncertainties into g bar and g obs as if they were random and uncorrelated (Lelli et al. 2017), or varied both the nuisance and RAR parameters galaxy-bygalaxy, effectively assuming a different RAR for each galaxy (Li et al. 2018;Chae et al. 2020bChae et al. , 2021Chae et al. , 2022)).Assuming an underlying universal form for the RAR, a superior inference constrains global RAR parameters along with the local galaxy properties.The main advantage of this is that it propagates the prior distributions of the galaxy parameters as systematic rather than statistical uncertainties, thus capturing the correlations across rotation curves that fluctuations in these parameters induce.For example, a higher (lower) Υ than expected in a particular galaxy causes a higher (lower) g bar across its rotation curve, yet modelling it as a statistical uncertainty implicitly assumes that a fluctuation in Υ could scatter g bar up at one point and down at the next.The full inference also captures the degeneracies between the RAR and galaxy parameters, which have a non-trivial impact on the relation through the shape of the galaxy priors.This is the analysis I perform here.
Although conceptually simple, the full inference is technically challenging because it implies a vastly higherdimensional parameter space than the simplified versions.The analysis of Lelli et al. (2017) has two parameters (a0, σint), while that of Li et al. (2018) has four (D, i, Υ disk , Υ bulge ) repeated N = 147 times for N galaxies.(Chae et al. 2020b additionally sample Υgas and eN.) Li et al. and Chae et al. cannot accommodate parameters that couple the galaxies, so fix a0 = 1.2 × 10 −10 m s −2 a priori and can at best reconstruct σint post-hoc from the distribution of residuals, thus neglecting its degeneracy with the other variables.The full inference has up to 6N +n+2 parameters (a0, σint, N ×eN, where n = 31 is the number of galaxies with bulges.Thus, although the total number of parameters that I sample is only slightly larger than Chae et al., the fact that I sample them together while Li et al. and Chae et al. split them by galaxy makes for a qualitatively different analysis, capable of mapping out the degeneracies between all parameters and inferring a0 and σint.915 parameters is indeed beyond many sampling methods, but routine for Hamiltonian Monte Carlo.This will enable a robust determination of the RAR parameters for arbitrary priors and assumptions about the underlying functional form.g bar and g obs transformed according to the best-fit galaxy parameter values (Fig. 1) reveals the RAR that underlies the sampling distributions of those parameters.
The structure of this paper is as follows.In Sec. 2 I describe the Spitzer Photometry and Accurate Rotation Curves (SPARC) data and selection criteria I employ.Sec. 3 gives the methodology, including the likelihood model, priors, treatment of the galaxy parameters and details of the sampler.The results are presented in Sec. 4. Sec. 5 discusses the broader ramifications of the study, remaining systematic uncertainties and useful further work, while Sec.6 concludes.Throughout, log has base 10 and accelerations are given in 10 −10 m s −2 unless otherwise stated.

OBSERVATIONAL DATA
I analyse the SPARC sample (Lelli et al. 2016),1 comprising 175 rotation curves from the literature with Spitzer photometry at 3.6µm.I apply the quality cuts recommended by Lelli et al. (2017), removing galaxies with quality flag 3 (indicating large asymmetries, non-circular motions and/or offsets between stellar and Hi distributions) or maximum a priori i < 30 deg, and points for which the quoted fractional uncertainty on the observed rotation velocity is greater than 10 per cent.This leaves 2696 points from 147 galaxies, of which all have mass in a stellar disk but only 31 have mass in a central bulge.
Distances are determined by a variety of methods with a corresponding range of uncertainties (Lelli et al. 2016), while the inclinations are estimated from tilted-ring fits to the velocity fields.I use these as Gaussian priors in the inference.The total luminosity at 3.6 µm, L3.6, is well-measured but its uncertainty is quoted so I include it as a Gaussian prior for completeness and to eliminate statistical uncertainty in the independent (g bar ) direction which complicates the inference (see Sec. 3.2).I follow the SPARC convention that L3.6 is calculated using the maximum a priori distance for each galaxy, D, and hence does not scale with D. Similarly, the uncertainty on L3.6, δL3.6, comes purely from the uncertainty on the flux and does not include a contribution from the distance uncertainty.The disk and bulge mass-tolight ratios, Υ disk and Υ bulge , are believed to be ∼ 0.5 and ∼ 0.7 respectively, with a ∼ 25 per cent uncertainty (Meidt et al. 2014;McGaugh & Schombert 2014;Lelli et al. 2016).I use these as lognormal priors, which are marginally favoured over Gaussian given the way the parameters are determined (S.McGaugh and F. Lelli, priv. comm.).δL3.6 is sufficiently small for it not to make a difference whether it is modelled as normal or lognormal.
With Hi mass measured, a correction factor must be applied to calculate the total gas mass and hence the gas contribution to g bar .The fiducial SPARC analysis uses a conversion factor of 1.33 (accounting for primordial helium), but a more accurate determination includes a scaling of the hydrogen fraction with the stellar mass M * of the galaxy (McGaugh et al. 2020): where X(M * ) = 0.75 − 38.2 M * /(1.5 × 10 24 M⊙) 0.22 . (2) As the Hi mass has already been scaled by 1.33 in SPARC, I define where overbar denotes maximum a priori value.(M * must be determined after sampling L3.6, Υ disk and Υ bulge .)This scales Υgas relative to the value assumed in SPARC when calculating Vgas, as Υ disk and Υ bulge do for V disk and V bul .The results are not significantly altered compared to Ῡgas = 1.
Υgas is given a lognormal prior with 10 per cent width (Lelli et al. 2016).g pred obs = g bar /2 + g 2 bar /4 + g bar a0. (6) Although in tension with Solar System measurements this function is highly successful for galaxy dynamics (Famaey & McGaugh 2012), and may readily be tweaked to circumvent local constraints without appreciably altering its larger-scale behaviour.One such modification is the "RAR IF" of Lelli et al. (2017), which I have checked yields almost identical results to the Simple IF.The IF currently has no physical significance and must be constrained empirically (Milgrom 2016;Famaey & McGaugh 2012).
The reason I use the Simple IF is that the second function I consider is designed to reduce to it in the zero-externalfield limit.This is the EFE formula for the nonrelativistic AQUAdratic Lagrangian (AQUAL; Bekenstein & Milgrom 1984) theory of MOND designed in Chae & Milgrom (2022): 2 To convert between normal and lognormal distributions I use the full equations relating their means and standard deviation where a tilde indicates the lognormal.The uncertainties are sufficiently small in most cases for this not to differ appreciably from the more common first-order approximation.
Figure 1.The underlying RAR of the SPARC sample (blue) is obtained by transforming g bar and g obs according to the best-fit galaxy parameters, in this case those at the median of the posterior for the inference with intrinsic scatter but without EFE.The Simple IF fit with a 0 in its 2σ allowed range is overplotted in red, and the median errorbar size, deriving solely from the statistical uncertainty in V obs , is shown as a magenta bar in the lower right.The standard "prior RAR", where the galaxy parameters take their maximum a priori values, is shown in faded grey.
, where eN ≡ gext/a0 describes the strength of the external field at the galaxy in question.The EFE arises in most formulations of MOND due to the theory's nonlinearity: the strong equivalence principle is violated because the acceleration of a system as a whole cannot be transformed away in calculation of its internal motions.This implies that otherwise identical galaxies in different gravitational environments have different kinematics.A stronger external field pushes the system towards the Newtonian regime by reducing the gravitational boost of MOND, causing a downturn in the RAR at low g bar where gext can be a non-negligible fraction of g obs .While several fitting formulae for the EFE exist (e.g.Banik & Zhao 2015;Haghi et al. 2019;Zonoozi et al. 2021), Eq. 7 is the most sophisticated in allowing for variable disk thickness and scale length-and the orientation of the field relative to the disk axis through azimuthal averaging-and has been shown to yield good agreement with the SPARC data (Chae et al. 2022;Chae 2022).It should be borne in mind however that this does not make it correct in general.I consider both the case of eN as a global parameter describing the average external field over the sample, and as a parameter varying galaxy-by-galaxy to describe their separate local environments.In the former case I use a uniform prior sufficiently broad to enclose the full posterior; in the latter, where there is insufficient information in the data for a meaningful constraint on eN, I impose a prior based on the environmental field estimates of the SPARC galaxies from Desmond et al. (2018);Chae et al. (2021).These are determined entirely independently of the SPARC data by summing contributions to the gravitational field from the baryonic masses of surrounding objects, including a sophisticated treatment of survey incompleteness and other missing mass.
As my fiducial analysis I use the results assuming that missing baryons are strongly clustered around visible objects ("maximum clustering") because this is expected in MOND and was shown in Chae et al. (2021Chae et al. ( , 2022) ) to give good agreement with the SPARC rotation curves.I also consider an "average clustering" model that assumes a prior distribution midway between the "max clustering" and "no clustering" (missing baryons uncorrelated with visible objects) results, with a width given by half the difference between the two.This systematic uncertainty is larger than the statistical uncertainty in either clustering case separately.The most precise calculation of Chae et al. ( 2021) uses data from the Sloan Digital Sky Survey and hence is only valid within the footprint of that survey, which includes 90 galaxies in my sample.For the remaining 57 I take ēN to be the median eN over all SPARC galaxies (in the corresponding clustering model), with an uncertainty twice the median uncertainty for all SPARC galaxies.This corresponds to a conservatively wide prior for galaxies without object-specific prior information, while still leveraging information on the eN distribution across the population.Combined with the no-EFE case (Eq.6), these EFE models ought roughly to span the space of possible EFE behaviour and hence indicate the level of systematic uncertainty that the unknown EFE behaviour induces.
Table 1 summarises the free parameters of the inference and their priors.

Inference procedure
The parameters inferred in the fiducial model are a0, σint, 147×eN, 147×D, 147×i, 147×L3.6,147×Υ disk , 147× Υgas and 31 × Υ bulge .At any point in parameter space I calculate g bar and g obs as where Vgas, V disk , and V bul are the velocities generated by the gas, disk and bulge, V obs is the observed velocity and r is the galactocentric radius.These are as quoted in the SPARC database, i.e. assuming D = D, i = ī, L3.6 = L3.6 and all Υ = 1.Vgas |Vgas| is used rather than V 2 gas in g bar to account for the possibility of central "holes" in the gas distribution which can cause the gravitational field sourced by the gas to point outwards.Note that g bar is independent of D because all of V 2 gas , V 2 disk , V 2 bul and r scale proportionally to D. The only remaining uncertainty to treat as statistical is the contribution of δV obs to g obs .3I assume this is lognormal, so that δ log(g obs ) = log(1 + (δg obs /g obs ) 2 )/ ln(10) (10) where δg obs /g obs = 2 δV obs /V obs .
I then use either Eq. 6 or 7 to calculate the predicted g obs at each g bar .σint simply adds in quadrature with δ log(g obs ), so the likelihood is: where ⃗ p is the parameter vector and j runs over the 2696 data points.This fiducial analysis assumes no statistical uncertainty on the velocities sourced by the gas, disk and bulge.However the calculation of V disk , V bul and Vgas in Lelli et al. (2016) made assumptions about the 3D geometry of these baryons, particularly in the thickness of the disk components.Variation may be expected to alter the baryon velocities at the ∼10-15 per cent level (F.Lelli, priv. comm.).I therefore also consider models in which these velocities are each assigned 10 per cent uncorrelated Gaussian uncertainties,4 which are propagated according to The application of these equations will be indicated by "boosted uncertainties".In this case, the presence of uncertainties in the x direction of the RAR plane introduces latent variables describing the true position of each point on the x-axis in the Bayesian hierarchical model.This makes the likelihood function for the parameters of interest alone ambiguous.Two approaches to remove the latent nuisance parameters without sampling them are to marginalise over them with a uniform prior, or to maximise the likelihood with respect to each of them (as a function of the other parameters in the inference) to produce a profile likelihood for the other parameters.These result in different maximum-likelihood points and parameter constraints.Tests on mock data (in agreement with literature results; Berger et al. 1999;Hadzhiyska et al. 2023) show that the marginalised likelihood recovers the correct intrinsic scatter and weakly biased shape parameters (e.g.a0 and eN) while the profile likelihood recovers unbiased shape parameters but can bias σint significantly low.As I am mainly interested in whether the boosted uncertainties allow for an intrinsic-scatter-free RAR I opt for the marginalised likelihood, which replaces Eq. 13 by δ log(g bar,j ) 2 (16) in Eq. 12 (for the derivation see e.g.sec.3.2 of Desmond et al. 2023).I perform the inference using the No U-Turns Sampler (NUTS; Hoffman & Gelman 2011) method of Hamiltonian Monte Carlo (HMC), as implemented in numpyro (Phan et al. 2019;Bingham et al. 2019).I initialise the sampler to the median of 20,000 points randomly drawn from the prior, which I find to yield good convergence behaviour.For each inference I concatenate 28 separate chains run in parallel, manually tuning the number of warmup and sampling steps to ensure that burn-in is complete and that there are enough effective samples for the Gelman-Rubin statistic (Gelman & Rubin 1992) to satisfy |r − 1| < 0.001.This requires ∼1000 warmup steps, ∼4000 sampling steps and takes ∼1 hour per model to run.

Validation with mock data
My analysis uses uniform priors on a0, σint and global eN, which are not reparametrisation invariant and cannot, pace popular opinion, be considered uninformative.Such priors are prone to contributing volume effects to the posterior, which can lead to significant biases when applied to parameters to which the likelihood is relatively insensitive (e.g.Hadzhiyska et al. 2023).In addition, the uncertainties and finite sample size lead to scatter in the maximum-likelihood parameters around the population values (sample variance).To assess the impact of these effects I analyse mock data generated by the following procedure: (1) Randomly sample the galaxy parameters from their prior distributions (2) Rescale the "observed" g bar according to Eq. 8 to calculate g true bar (3) Use either the no-EFE, global-EFE or max-clustering-EFE model, with some true a0, eN, σint, to calculate g true obs from the g true bar (4) Transform to observed g obs through Eq. 9, and hence to V obs assuming the same r values in the mock data as in the SPARC data (5) Replacing the SPARC V obs by these values, calculate the maximum-likelihood values of all parameters and run the inference to compute their posteriors (6) Repeat twice with different random seeds.
Choosing a0 = 1.2, σint = 0.05 dex and eN = 0.01 in the case with global EFE, the results are shown in Table 2.I find neither the maximum-likelihood parameters nor their posteriors to be significantly different to their true values, showing the above effects not to be important for the data and models under consideration.The eN values and galaxy parameters (not shown) are similarly unbiased.This gives confidence to proceed with the analysis of the real data.

RESULTS
Table 3 shows the median and 1σ uncertainty of the RAR parameters for each of the models considered, along with their maximum log-likelihood (ln( L)), maximum logposterior (ln( P )) and Bayesian information criterion (BIC) relative to the first model.The BIC should be taken as a very rough estimator only for the Bayesian evidence, both because the number of data points does not greatly exceed the number of parameters and because the parameter priors are not necessarily slowly varying at the maximum a-posteriori point.For example, replacing 2 ln( L) by 2 ln( P ) in the BIC formula changes ∆BIC to 1.44 for the "No scatter, global EFE" model, turning "decisive" evidence on the Jeffreys scale against the inclusion of eN to "barely worth mentioning."Regardless of this, a clear result from the goodness-of-fit statistics is that either non-zero intrinsic scatter or boosted uncertainties is strongly preferred, mainly through a large increase in the likelihood.The fiducial model, which I considered a priori to be most likely, is "Scatter, max-clustering EFE".
The value of a0 is slightly reduced if the RAR is assumed to possess intrinsic scatter, and more significantly increased by including the EFE with external field strength priors from the baryonic large-scale structure.This is because these priors somewhat increase eN relative to likelihood alone and there is a positive degeneracy between eN and a0.Differences in a0 between the models are up to a few times larger than their statistical uncertainties, and a naive averaging over all the models implies a0 = 1.19 ± 0.04 (stat) ± 0.09 (sys).
The intrinsic scatter is similar in all models in which it is included, except in the "boosted uncertainties" case.The results suggest that σint = 0.034 ± 0.001 (stat) ± 0.001 (sys) dex, with the significant caveat that this assumes the SPARC error model is reliable.A 10 per cent uncertainty on V disk , V bul and V disk -as may be expected from deviations from the assumed 3D baryon geometry-is sufficient to set the preferred σint to 0. The intrinsic scatter is not driven by the most egregious outliers: removing the 11 points with g bar > g obs and g obs < 2 after transformation for the model without EFE (see Fig. 1) reduces it only to 0.031 dex.
There is weak evidence for ⟨eN⟩ > 0 when it is given a uniform prior, as evidenced both by the inferred eN being consistent with 0 within ∼3σ and by the small gain in maximum posterior value across the chain (∆ ln( P )).Adding this parameter is not favoured by the BIC.The model perhaps most similar to the MOND expectation with global eN (not shown in Table 3) has σint = 0 and boosted uncertainties: in this case ⟨eN⟩ = 0.0024 ± 0.0009, again a 3σ "detection" but a slightly larger value.This model has ∆BIC= −3860, very similar to the case with σint and boosted uncertainties but without EFE, showing that again the addition of a global eN is not favoured.
The average prior values of eN over all the galaxies for the maximum-clustering and average-clustering cases are 0.0050 and 0.0018 respectively, while the posteriors average to 0.0050 and 0.0022 without σint and 0.0049 and 0.0020 with.This shows that when modelled galaxy-by-galaxy the data does not disfavour significant values of eN in agreement with the large-scale structure expectation, and in fact in both cases the maximum-likelihood value is increased over the case of no or global eN, significantly without intrinsic scatter and moderately with it.Prior evidence is however required to support the presence of the EFE in a model comparison sense as the BIC strongly disfavours the addition of 147 eN parameters.While the average-clustering EFE model gives a slightly higher L when including σint, it gives a lower L without it and a lower P in both cases.This constitutes weak evidence in favour of the maximum-clustering prior, but the data is far from sufficient to distinguish robustly between them.Fig. 1 uses Eqs. 8 and 9 to transform g bar and g obs according to the parameters at the median of the posterior for the inference with σint but without EFE (the model preferred by the BIC), including the 2σ model prediction.This illustrates the extreme tightness of the underlying (posterior) RAR, including relative to the traditional "prior RAR" shown in grey.The few outlying points are also outliers of the prior RAR, and are brought slightly closer to the line by the transformation.Note that the systematic uncertainty on the blue points due to their dependence on the galaxy nuisance parameters is not shown.
Fig. 2 shows an excerpt from the corner plot of the inference that also includes the EFE with a global eN.The posteriors on three parameters for galaxy 1 (D512-2) are compared to their maximum prior values (blue lines), with which they agree well.a0, eN and σint are not strongly degenerate with other parameters, while the full parameter space exhibits the degeneracies expected from Eqs. 8 and 9.The median distance to the SPARC galaxies is reduced by 1.7 per cent going from the inference with σint but without EFE to that with local eN and maximum-clustering prior, although within the latter inference D is positively correlated with eN on average within the chain.This is because at higher eN, g pred obs is lower at low g bar and hence closer to g obs at higher D (or i).
To show more generally the differences between the prior and posterior values of the galaxy parameters, Fig. 3 shows the distribution of normalised residuals for the model with intrinsic scatter and galaxy-by-galaxy eN with maximumclustering prior.For galaxy parameter X, the normalised residual is defined as where subscript "post" and "prior" denote the posterior and prior distributions, and "med" and "std" stand for median and standard deviation.Larger values of the residual indicate a posterior significantly shifted from the prior due to the influence of the likelihood.The width of the distributions therefore reflect the sensitivity of the data to the parameters: those with distributions sharply peaked at 0 such as L3.6 and Υgas are relatively unimportant so that med(Xpost) ≈ X, while those with very broad distributions have likelihood often peaked far from the prior centres.This is especially pronounced for Υ bulge and Υ disk which determine g bar in the galaxies' central regions.Υ bulge in particular has a slight negative offset to reduce the number of inner points for which g bar > g obs .The black dashed lines are what one would expect if Xpost scatters around X with a Gaussian distribution of width given by std(Xpost).Finally, I show in Fig. 4 smoothed distributions of the fractional uncertainties in the galaxy-specific parameters in the prior and posterior for the inference with intrinsic scat- Table 3. Constraints on RAR parameters and goodness-of-fit statistics for the models considered.For the maximum-clustering and average-clustering EFE models the quoted e N constraints describe the stacked posteriors over all galaxies.The final three columns are the maximum log-likelihood, maximum log-posterior and Bayesian information criterion relative to the first model.
ter and local eN with maximum-clustering prior.In most cases the RAR constraint has increased the precision with which the parameters are known.This is especially marked for the distance, where galaxies in the second mode of the prior distribution (those at ≲60 Mpc with only redshift distances; Lelli et al. 2016) are brought into a single posterior mode at ∼10 per cent uncertainty.This illustrates the utility of the RAR as a direct (i.e.redshift-independent) distance measurement method.There is analogous behaviour for eN, where the higher prior mode corresponds to galaxies outside the SDSS footprint which are assigned higher uncertainties (see Sec. 3.1); these become better known on applying the RAR constraint.The uncertainties on Υ disk and Υ bulge also fall markedly, partly to reduce the number of inner points with posterior probability at g bar > g obs . 5The analogues of Figs.1-4 for the other models are qualitatively similar, in line with the variation in their results shown in Table 3. Supplementary tables available at https://zenodo.org/record/7752545 (Desmond 2023) contain the mean, median and 1, 2 and 3σ constraints on the parameters for each of the models considered.The format is illustrated in Table 4, which shows the first five and last two rows for the model with σint and global eN.An example use of these is to resolve distance variations among galaxies in the Ursa Major cluster, all of which have D = 18 Mpc.A reference table contains the galaxy names at each index as well as the bulge luminosity and means and standard deviations of their Gaussian priors on D, i and L3.6, transcribed from the SPARC database.Plots of the prior and posterior constraints on D, i, L3.6, eN, Υ disk , Υ bulge and Υgas for all galaxies under each model are also included.

DISCUSSION
While the value of a0 that I infer (1.1 ≲ a0 ≲ 1.3) is in full agreement with literature results, the RAR intrinsic scatter σint ≈ 0.034 dex is significantly smaller.Lelli et al. (2017) quote a total scatter of 0.13 dex and argue that most of this comes from observational uncertainties; using their model and Eqs. 12 and 16 I find σint = 0.082 ± 0.003 dex.Li   (2018) quote an intrinsic scatter of 0.057 dex from residuals around their best-fit relation with fixed a0 = 1.2.A reduction in intrinsic scatter when marginalising over galaxy parameters is not guaranteed: although the denominator of Eq. 12 favours smaller σint, this is offset by lower prior probabilities of the galaxy parameters at values that bring the points closer to the theoretical line to reduce the exponent in Eq. 12.
The statistical uncertainties on g bar and g obs also fall greatly (the former to 0) when the galaxy parameters are inferred, so a lower total scatter of the points around the line need not translate into lower σint.That the preferred intrinsic scatter in the full analysis is extremely small hints towards the RAR being at base a practically monotonic correlation.
The question of the fundamentality of the RAR has im- portant ramifications for the mass discrepancy problem on galaxy scales.Lelli et al. (2017) claimed the relation to be law-like, a result supported here.It will be challenging for simulations or semi-analytic models in ΛCDM to reproduce the tightness of the underlying RAR given that even the scatter of the prior RAR, ∼0.1 dex, is non-trivial (Di Cintio & Lelli 2016;Desmond 2017;Keller & Wadsley 2017;Ludlow et al. 2017).It would be interesting to explore this explicitly using modern high-resolution cosmological hydrodynamical ΛCDM simulations such as TNG-50 (Pillepich et al. 2019;Nelson et al. 2019) and NewHorizon (Dubois et al. 2021), which would provide a stringent test of those models.
It is important to bear in mind that SPARC is only a small fraction of the data pertinent to the RAR: one may also use ultra-diffuse galaxies (e.g.Freundlich et al. 2022), local dwarf spheroidals (e.g. McGaugh & Wolf 2010;Mc-Gaugh & Milgrom 2013), early-type galaxies (e.g.Rong et al. 2018;Chae et al. 2020a), low-acceleration regions including the outer Milky Way (Oman et al. 2020), stacked weak lens-ing (Brouwer et al. 2021) and groups and clusters of galaxies (Chae et al. 2019(Chae et al. , 2020a;;Chan & Del Popolo 2020;Tian et al. 2020;Pradyumna & Desai 2021;Gopika & Desai 2021).Some of this data appears to deviate from the MOND expectation.SPARC is however one of the few datasets with uncertainties under sufficient control for the present analysis to be feasible and meaningful.Future work may extend it to new regimes.
Naïvely one expects σint = 0 in MOND, but in fact this is only true in modified inertia formulations in the limit of perfectly circular orbits.Even in such models a scatter is introduced by deviations from circularity because the dynamics of an object depends on its entire past trajectory (Milgrom 2011;Milgrom 2022).This may already be sufficient to account for the 0.034 dex (8 per cent) scatter present in the underlying RAR, suggesting that modified inertia, which predicts most directly the algebraic relation that I fit, is viable.Additional scatter is present in modified gravity formulations where the algebraic MOND relation holds only in spherical symmetry (Famaey & McGaugh 2012), a condition clearly violated in disk galaxies.It would be interesting to quantify the scatter introduced by these effects in SPARC-like galaxies, further testing the MOND paradigm and providing a novel way to distinguish between the modified gravity and modified inertia interpretations (Petersen & Lelli 2020;Chae 2022).The model of constant mass-to-light for the disk and bulge may also be overly simplistic, with radial dependence parameters able to soak up some of the remaining scatter.
The presence or absence of the EFE is controversial within the MOND literature, with some studies claiming strong evidence for it (McGaugh & Milgrom 2013;Haghi et al. 2019;Chae et al. 2020b) and others strong evidence against (Hernandez et al. 2021;Hernandez & Lara-D I 2019;Freundlich et al. 2022).My work does not resolve this issue: there is weak evidence for a positive average external field strength across the sample, while inferring it galaxy-by-galaxy with a prior from independent measurements of environment improves the likelihood but is not favoured by the Bayesian information criterion.An important caveat is that the fitting formula I use (Eq.7) was only designed for the outer regions of rotation curves (Chae & Milgrom 2022)-while I am applying it to them in their entirety-and is only valid within the AQUAL model.Different EFE formulae may give different constraints on eN, and hence a0 with which it is degenerate, but would not be expected to alter σint appreciably.(A scatter in the effect of the EFE, e.g.due to variably internal and external fields, may however reduce the true intrinsic scatter.)The environmental priors are also highly uncertain due to the possibility of clustered unseen baryonic mass.Further work on modelling the EFE and constraining the external field is therefore required to reach a definitive conclusion concerning the existence of the EFE in galaxy dynamics.As MOND is currently an effective model only, it may be that the underlying theory gives a mass or scale dependence to the EFE which can reconcile seemingly discrepant results.
Besides calibrating the RAR my work provides strong constraints on the properties of the SPARC galaxies under the assumption that the underlying RAR is as I model it.Summaries of these constraints are made public in online tables to facilitate future studies using SPARC.One such application is to use the RAR as a direct distance probe.The current study calibrates the relation using a sample with informative distance priors.The distance to any galaxy with a (partially) resolved rotation curve may be inferred by fitting it to the calibrated RAR, including marginalisation over the other relevant properties of the galaxy but not necessarily the parameters of the RAR itself.This is analogous to the wellestablished Tully-Fisher and Fundamental Plane methods, but achieves higher precision (∼10 per cent uncertainty on D rather than 20-25 per cent; Tully et al. 2023) at the cost of requiring resolved kinematics.This may readily be achieved for a large sample of galaxies using the high spatial resolution of upcoming instruments such as the Square Kilometer Array.Note however that the RAR may evolve with redshift, e.g.due to time-dependent a0 in MOND or evolution of galaxy and halo density profiles in ΛCDM (Keller & Wadsley 2017;Paranjape & Sheth 2021), which would necessitate recalibration of the relation when this effect kicks in.My analysis also supplies enhanced kinematic inclinations, as well as constraints on mass-to-light ratios which may be correlated with other galaxy properties to advance understanding of stellar populations and galaxies' gas content.
Outliers of the underlying RAR may either be individual rotation curve points with large residuals from the best-fit line (e.g. in Fig. 1) or entire galaxies with parameters strongly shifted from their prior centres to achieve a good fit (readily visible in the supplementary figures).Studying these on a case-by-case basis may help identify peculiar galactic features and signpost the need for more sophisticated modelling.For example, NGC 2915 seems to require very high Υ disk ≈ 1.2 in all models.This is a starburst dwarf galaxy with significant radial motion in the inner regions, the modelling of which affects the entire rotation curve.It has a highly complex structure and is likely not in dynamical equilibrium towards the centre, while further out there is a strong warp (Meurer et al. 1996;Elson et al. 2010Elson et al. , 2011a,b;,b;Tang et al. 2022).The ability to correlate such properties with the results of this analysis across the sample would lend weight to the interpretation of the underlying RAR as fundamental.
Provided the fitting functions are good, the constraints on galaxy parameters do not assume MOND any more than modelling the Tully-Fisher relation as a power-law.Both may be thought of as empirical descriptions of the data without consideration of their theoretical significance.Recently, however, Desmond et al. (2023) have challenged the optimality of MOND functions (those with Newtonian and deep-MOND or EFE-driven regimes) for fitting the RAR, finding that the majority of functions that most efficiently compress the SPARC data tend to constant g obs at low g bar .These functions and their parameters have no known theoretical significance, but would lead to different "underlying" relations and hence galaxy parameter constraints.This systematic uncertainty in e.g.distance measurements could be assessed by repeating the present inference with one or more of these functions.It is however not known how the results of Desmond et al. would be affected by marginalising over the galaxy parameters separately for each function as done here, which would alter the function ranking, or by extending Exhaustive Symbolic Regression (Bartlett et al. 2022) to higher complexity.
I see three technical ways in which this inference could be improved.The first is to separate δV obs into a truly statistical and systematic part, including its full covariance across the galaxies' rotation curves.The average δ log(V obs ) across the sample is 0.03 dex, so more accurate modelling of this has the potential to alter the best-fit value of σint, which is similar.In particular, an important contribution to δV obs is from the difference in velocity between the approaching and receding sides of the disk, which is likely strongly correlated over r.This may itself be correlated with inclination, which in detail may vary across the disk.The second is to characterise better the molecular gas and 3D baryon geometry.The former could be done by estimating (or sampling) MH 2 from the M * − MH 2 relation of a population of similar galaxies and adding a contribution to V bar from a thin H2 disk, and the latter by solving Poisson's equation for different assumptions about disk thickness, bulge oblateness and potential asymmetries.For the SPARC galaxies these uncertainties are however highly subdominant to those that I model explicitly.The third is to use a Jeffreys rather than uniform prior for a0, eN and σint, thus eliminating potential volume effects.The mock tests of Sec.3.3 show these not to bias the results significantly.A Jeffreys prior would however enable the inference of eN galaxy-by-galaxy without importing large-scale structure information; when I tried this on mock data using a uniform prior I found a0 to be biased high, presumably due to the allowed and poorly-constrained volume towards high eN.Another extension would be to try a more complex intrinsic scatter model; Li et al. (2018) for example find a superposition of two Gaussians to fit the residuals better than one, perhaps reflecting the separate formal error and kinematic asymmetry contributions to V obs .

CONCLUSION
I have uncovered underlying RARs in the SPARC data by fitting the parameters of the Simple interpolating function, with and without intrinsic scatter and the external field effect, simultaneously with all relevant galaxy properties.The preferred intrinsic scatter is very small, 0.034±0.002dex, and additional plausible uncertainties are capable of reducing this to 0. The acceleration constant is in the range 1.1 ≲ a0 ≲ 1.3, in good agreement with literature results.I find weak evidence for the external field effect using an average external field strength over the full sample with a uniform prior.Allowing the field strength to vary galaxy-by-galaxy with a prior from large-scale structure observations improves the overall likelihood of the data but is not favoured by the Bayesian information criterion.My results suggest near-monotinicity and a high degree of regularity in the the RAR, providing a fresh challenge to galaxy formation models.The constraints I produce on the SPARC galaxies' parameters-distance, inclination, luminosity and disk, bulge and gas mass-to-light ratios-are the most precise to date (although subject to systematic error if the underlying RAR is not as I model it).I publicly release summaries of all posteriors for analyses that may benefit from this information.

DATA AVAILABILITY
The mean, median and 1, 2 and 3σ confidence intervals of the parameters for all models are available at https: //zenodo.org/record/7752545.The remaining data generated here, including the full HMC chains to explore degeneracies, will be made available on reasonable request.
functions to the RAR.The first is the "Simple interpolating function (IF)" (Famaey & Binney 2005):

Figure 2 .
Figure 2. Corner plot of selected parameters from the inference including intrinsic scatter and global e N .For the distance (in Mpc), inclination (in degrees) and disk mass-to-light ratio of the first galaxy (D512-2), the maximum a priori values are shown by the blue lines.

Figure 3 .
Figure3.The distribution of normalised residuals (Eq.17) of the galaxy-specific parameters in the inference including intrinsic scatter and local e N with maximum-clustering prior.The mirroring of some of the curves about the x-axis is for visual clarity only.A standard normal distribution-corresponding to Xpost scattering around X as a Gaussian with width given by std(Xpost)-is shown in dashed black.

Figure 4 .
Figure 4. Distribution of prior and posterior fractional uncertainties on galaxy-specific parameters for the inference with scatter and local e N with maximum-clustering prior.In the lower panels the vertical dashed lines are the fractional prior uncertainties (25 and 10 per cent), which are the same for all galaxies.

Table 1 .
The free parameters of the inference and their priors.N t denotes a truncated normal with lower and upper limits given by the final two arguments.

Table 4 .
et al.Excerpt of posterior summaries for the model with intrinsic scatter and global e N .The full tables for all models are available online.