# On the functional form of the radial acceleration relation

Harry Desmond, Deaglan J. Bartlett, Pedro G. Ferreira

► To cite this version:

Harry Desmond, Deaglan J. Bartlett, Pedro G. Ferreira. On the functional form of the radial acceleration relation. Monthly Notices of the Royal Astronomical Society, 2023, 521 (2), pp.1817-1831. 10.1093/mnras/stad597. hal-03960038

HAL Id: hal-03960038

https://hal.science/hal-03960038

Submitted on 3 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the functional form of the radial acceleration relation

Harry Desmond<sup>1,1*</sup> Deaglan J. Bartlett<sup>2,3*</sup> and Pedro G. Ferreira<sup>3</sup>

Accepted 2023 February 21. Received 2023 February 20; in original form 2023 January 12

## ABSTRACT

We apply a new method for learning equations from data – *Exhaustive Symbolic Regression* (ESR) – to late-type galaxy dynamics as encapsulated in the radial acceleration relation (RAR). Relating the centripetal acceleration due to baryons, $g_{bar}$, to the total dynamical acceleration, $g_{obs}$, the RAR has been claimed to manifest a new law of nature due to its regularity and tightness in agreement with Modified Newtonian Dynamics (MOND). Fits to this relation have been restricted by prior expectations to particular functional forms, while ESR affords an exhaustive and nearly prior-free search through functional parameter space to identify the equations optimally trading accuracy with simplicity. Working with the SPARC data, we find the best functions typically satisfy $g_{obs} \propto g_{bar}$ at high $g_{bar}$, although the coefficient of proportionality is not clearly unity and the deep-MOND limit $g_{obs} \propto \sqrt{g_{bar}}$ as $g_{bar} \to 0$ is little evident at all. By generating mock data according to MOND with or without the external field effect, we find that symbolic regression would not be expected to identify the generating function or reconstruct successfully the asymptotic slopes. We conclude that the limited dynamical range and significant uncertainties of the SPARC RAR preclude a definitive statement of its functional form, and hence that this data alone can neither demonstrate nor rule out law-like gravitational behaviour.

**Key words:** methods: data analysis – galaxies: kinematics and dynamics – dark matter.

## 1 INTRODUCTION

Kinematic measurements of galaxies relate their visible and dynamical masses, affording constraints on the distribution of dark matter and/or the behaviour of gravity. These measurements are simplest to perform for late-type galaxies supported predominantly by rotation, as the enclosed dynamical mass may be calculated from the centripetal acceleration and the law of gravity. Such studies have revealed a striking correlation between the enclosed baryonic and total dynamical mass assuming Newtonian gravity, dubbed the mass discrepancy–acceleration (Sanders 1990; McGaugh 2004) or radial acceleration relation (RAR; Lelli et al. 2017). It has been claimed that the RAR indicates that at high accelerations the Newtonian dynamical mass follows the baryonic mass (indicating little dark matter and the validity of Newtonian mechanics), while as acceleration drops below a new constant of nature $g_0 \approx 10^{-10} \text{ m s}^{-2}$, the dynamical mass increasingly exceeds the baryonic mass in a regular way.

One may attempt to understand these observations from either a dark matter or modified gravity perspective. In $\Lambda$CDM, the difference between the dynamical and baryonic mass is due to the dark matter that makes up most of the mass of the galaxy. The RAR must therefore be explained by the relative distributions of dark and visible mass established by the process of galaxy formation. Interactionless cold

dark matter is influenced only gravitationally by the baryonic mass so the emergence of the RAR must be somewhat fortuitous; it is not established directly by a baryon–dark matter coupling (although see Blanchet & Le Tiec 2008; Berezhiani & Khoury 2015; Famaey, Khoury & Penco 2018 for alternative ideas). In contrast, the modified gravity (or modified inertia) interpretation posits a breakdown of Newtonian mechanics at low acceleration so that the dynamical mass inferred by a Newtonian analysis is not the true dynamical mass of the galaxy. The prototypical instantiation of this idea is Modified Newtonian Dynamics (MOND; Milgrom 1983a, c, b), in which the kinematic acceleration $g_{obs}$ follows the square root of the Newtonian acceleration $g_{bar}$ in the weak-field regime. This enables the total dynamical mass of the galaxy to remain equal to the baryonic mass across galaxies’ rotation curves, eliminating the need for dark matter in them. The MOND paradigm attempts to dispense with dark matter entirely, and has cosmologically viable relativistic extensions (most recently Skordis & Złośnik 2021). It is reviewed in Famaey & McGaugh (2012) and Banik & Zhao (2022).

Central to the dark matter–modified gravity debate in the context of galaxy dynamics is the functional form of the RAR. This is because MOND makes a very specific prediction (absent the external field effect: $g_{obs} = g_{bar}$ in the high-acceleration ‘Newtonian regime’ and $g_{obs} \propto g_{bar}^{1/2}$ in the low-acceleration ‘deep-MOND regime’) while dark matter could accommodate a range of possibilities depending on the effect of galaxy formation on halo density profiles, which remains highly uncertain (e.g. Duffy et al. 2010; Macciò et al. 2012; Keller & Wadsley 2017; Ludlow et al. 2017; Navarro et al. 2017; Tenneti et al. 2018; Grudić et al. 2020). The only potentially

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

<footnote n="2"><sup>2</sup>CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>* E-mail: harry.desmond@physics.ox.ac.uk (HD); deaglan.bartlett@physics.ox.ac.uk (DJB)</footnote>
# On the functional form of the radial acceleration relation

Harry Desmond<sup>1,1*</sup> Deaglan J. Bartlett<sup>2,3*</sup> and Pedro G. Ferreira<sup>3</sup>

Accepted 2023 February 21. Received 2023 February 20; in original form 2023 January 12

## ABSTRACT

We apply a new method for learning equations from data – *Exhaustive Symbolic Regression* (ESR) – to late-type galaxy dynamics as encapsulated in the radial acceleration relation (RAR). Relating the centripetal acceleration due to baryons, $g_{bar}$, to the total dynamical acceleration, $g_{obs}$, the RAR has been claimed to manifest a new law of nature due to its regularity and tightness in agreement with Modified Newtonian Dynamics (MOND). Fits to this relation have been restricted by prior expectations to particular functional forms, while ESR affords an exhaustive and nearly prior-free search through functional parameter space to identify the equations optimally trading accuracy with simplicity. Working with the SPARC data, we find the best functions typically satisfy $g_{obs} \propto g_{bar}$ at high $g_{bar}$, although the coefficient of proportionality is not clearly unity and the deep-MOND limit $g_{obs} \propto \sqrt{g_{bar}}$ as $g_{bar} \to 0$ is little evident at all. By generating mock data according to MOND with or without the external field effect, we find that symbolic regression would not be expected to identify the generating function or reconstruct successfully the asymptotic slopes. We conclude that the limited dynamical range and significant uncertainties of the SPARC RAR preclude a definitive statement of its functional form, and hence that this data alone can neither demonstrate nor rule out law-like gravitational behaviour.

**Key words:** methods: data analysis – galaxies: kinematics and dynamics – dark matter.

## 1 INTRODUCTION

Kinematic measurements of galaxies relate their visible and dynamical masses, affording constraints on the distribution of dark matter and/or the behaviour of gravity. These measurements are simplest to perform for late-type galaxies supported predominantly by rotation, as the enclosed dynamical mass may be calculated from the centripetal acceleration and the law of gravity. Such studies have revealed a striking correlation between the enclosed baryonic and total dynamical mass assuming Newtonian gravity, dubbed the mass discrepancy–acceleration (Sanders 1990; McGaugh 2004) or radial acceleration relation (RAR; Lelli et al. 2017). It has been claimed that the RAR indicates that at high accelerations the Newtonian dynamical mass follows the baryonic mass (indicating little dark matter and the validity of Newtonian mechanics), while as acceleration drops below a new constant of nature $g_0 \approx 10^{-10} \text{ m s}^{-2}$, the dynamical mass increasingly exceeds the baryonic mass in a regular way.

One may attempt to understand these observations from either a dark matter or modified gravity perspective. In $\Lambda$CDM, the difference between the dynamical and baryonic mass is due to the dark matter that makes up most of the mass of the galaxy. The RAR must therefore be explained by the relative distributions of dark and visible mass established by the process of galaxy formation. Interactionless cold

dark matter is influenced only gravitationally by the baryonic mass so the emergence of the RAR must be somewhat fortuitous; it is not established directly by a baryon–dark matter coupling (although see Blanchet & Le Tiec 2008; Berezhiani & Khoury 2015; Famaey, Khoury & Penco 2018 for alternative ideas). In contrast, the modified gravity (or modified inertia) interpretation posits a breakdown of Newtonian mechanics at low acceleration so that the dynamical mass inferred by a Newtonian analysis is not the true dynamical mass of the galaxy. The prototypical instantiation of this idea is Modified Newtonian Dynamics (MOND; Milgrom 1983a, c, b), in which the kinematic acceleration $g_{obs}$ follows the square root of the Newtonian acceleration $g_{bar}$ in the weak-field regime. This enables the total dynamical mass of the galaxy to remain equal to the baryonic mass across galaxies’ rotation curves, eliminating the need for dark matter in them. The MOND paradigm attempts to dispense with dark matter entirely, and has cosmologically viable relativistic extensions (most recently Skordis & Złosnik 2021). It is reviewed in Famaey & McGaugh (2012) and Banik & Zhao (2022).

Central to the dark matter–modified gravity debate in the context of galaxy dynamics is the functional form of the RAR. This is because MOND makes a very specific prediction (absent the external field effect: $g_{obs} = g_{bar}$ in the high-acceleration ‘Newtonian regime’ and $g_{obs} \propto g_{bar}^{1/2}$ in the low-acceleration ‘deep-MOND regime’) while dark matter could accommodate a range of possibilities depending on the effect of galaxy formation on halo density profiles, which remains highly uncertain (e.g. Duffy et al. 2010; Macciò et al. 2012; Keller & Wadsley 2017; Ludlow et al. 2017; Navarro et al. 2017; Tenneti et al. 2018; Grudić et al. 2020). The only potentially

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

<footnote n="2"><sup>2</sup>CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Astrophysics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote>* E-mail: harry.desmond@physics.ox.ac.uk (HD); deaglan.bartlett@physics.ox.ac.uk (DJB)</footnote>
unambiguous prediction is that the RAR tends to $g_{\text{obs}} = \Omega_m/\Omega_b g_{\text{bar}}$ at radii sufficiently large to encompass the cosmic baryon fraction, but it is unclear where or even if this occurs in galaxies. Thus, while the $\Lambda$CDM prediction for the full RAR can be tested only by applying potentially restrictive priors on galaxy formation effects (Di Cintio & Lelli 2016; Desmond 2017; Paranjape & Sheth 2021), a more direct route towards informing the dark matter-modified gravity debate is to test the MOND prediction, specifically the limiting behaviour at $g \ll g_0$ and $g \gg g_0$, the small intrinsic scatter and the lack of residual correlations.

Here, we focus on the asymptotic behaviour. This can be assessed to some extent by fitting a functional form with free power-law slopes at both ends (Lelli et al. 2017), but this assumes that the slope tends to a constant at each end and restricts to a specific part of the functional parameter space for which this is the case. These are in question when assessing the accuracy of the MOND prescription. A fully satisfactory fit should therefore make no such assumptions, eliminating potential confirmation bias and testing without any priors the assertion that the RAR implies no dynamically relevant dark matter at high $g$ and the deep-MOND limit at low $g$. We accomplish this here by means of a novel regression algorithm dubbed *Exhaustive Symbolic Regression* (ESR; Bartlett, Desmond & Ferreira 2022b), and hence assess the degree to which the RAR supports the tenets of MOND. Within the MOND paradigm, this method also enables optimization of the ‘interpolating function’ (IF) $g_{\text{obs}} = F(g_{\text{bar}})$ between the two stipulated limits.

The structure of the paper is as follows. In Section 2, we describe the RAR data that we use, and in Section 3, our algorithm for generating functions and assessing their aptitude for describing the data. Section 4 presents the results. In Section 5, we discuss the broader ramifications, potential remaining uncertainties and ways in which the programme could be furthered in the future. Section 6 concludes; full details on ESR are given in the companion paper Bartlett et al. (2022b). Units not explicitly given are $10^{-10} \text{ m s}^{-2}$, and all logarithms are natural.

## 2 OBSERVATIONAL DATA

We use the SPARC data set (Lelli, McGaugh & Schombert 2016),<fnref n="1" /> a compilation of 175 rotation curves from the literature combined with Spitzer 3.6 µm photometry. We apply the same quality cuts as the RAR study of Lelli et al. (2017), removing galaxies with quality flag 3 (indicating large asymmetries, non-circular motions and/or offsets between stellar and HI distributions), and those with inclinations $i < 30^{\circ}$ 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.

## 3 METHOD

We describe our method for generating and assessing trial functions in Section 3.1, and our likelihood function in Section 3.2. In Section 3.3, we outline our criteria for assessing whether a function displays MOND-like behaviour.

### 3.1 Exhaustive symbolic regression

While algorithms for symbolic regression (SR) – the search for good functional descriptions of a data set – are becoming mature, they

remain fallible (La Cava et al. 2021). Unless the generating function of the data is known at the outset (in which case SR is not required), it is not possible to determine whether any SR algorithm has uncovered the best function. This motivated us to develop ESR which, given a set of basis functions, produces and evaluates every possible function up to a given complexity of equation, defined here as the number of nodes in its tree representation. This enables a brute-force solution to relatively simple problems and provides a touchstone for assessing the results of stochastic algorithms at higher complexity. As shown in detail in Bartlett et al. (2022b), stochastic searches regularly fail to the best functions at even moderate complexity ~6, so that we would not be confidant of obtaining the functional form of the RAR through any algorithm besides ESR.

Presented in full in the companion paper, ESR has two main steps: i) generating, and optimising the parameters of, all functions up to a given complexity, and ii) ranking these functions using an information-theoretic metric combining accuracy and simplicity. For part i, the steps are:

(i) Generate all possible trees containing a given number of nodes (equal to the complexity of functions considered).

(ii) Generate the complete set of such functions by decorating these trees with all permutations of the operators from the operator list specified in advance, utilising the constraints on the arity of the operator that can occupy a given node.

(iii) Simplify the functions and remove duplicates. Variants of the same function (e.g. $x(x + \theta_0)$ and $x^2 + \theta_0 x$) are however retained as these may have different model complexities (used in step ii, below). For each unique function the variant is retained that minimizes this.

(iv) Determine the values of the free parameters appearing in the functions that maximize the likelihood of the data (see Section 3.2).

(v) Repeat for all complexities under consideration.

The only degrees of freedom in this procedure are the maximum complexity considered (here set at 9 as higher complexity is computationally prohibitive) and the set of operators of which the functions are composed. Here we choose:<fnref n="2" />

(i) **Nullary**: $g_{\text{bar}}, \theta$

(ii) **Unary**: exp, sqrt, square, inv

(iii) **Binary**: $+$, $-$, $*$, $/$, pow

where $\theta$ is a free parameter. We implicitly take the absolute value of the argument of any square root or power.

The result of this procedure is a list of all functions up to the maximum complexity (of which there are $2.24 \times 10^7$), along with the parameter values that maximize the likelihood of the RAR data. As in regular regression, using the maximum likelihood as the model selection criterion would favour overfitting, whereby a function fits the data near-perfectly but generalizes or extrapolates poorly. To remedy this, SR typically uses two-objective optimization, where the second objective is the ‘simplicity’ of the function. In the absence of a metric for trading accuracy (the first objective) with complexity, optimal functions form a ‘Pareto front’ where accuracy cannot be increased without reducing simplicity and vice versa. Simplicity has been defined analogously to model complexity (the number of nodes in the tree representation; e.g. in PySR; Cranmer et al. 2020), among others, but such definitions are typically

<footnote n="1">¹http://astroweb.cwru.edu/SPARC/</footnote>

<footnote n="2">²Many of these operators can manifestly be constructed from combinations of others. We include these repetitions to simplify the function trees, allowing a greater range of expressions up to the maximum complexity.</footnote>
unambiguous prediction is that the RAR tends to $g_{\text{obs}} = \Omega_m/\Omega_b g_{\text{bar}}$ at radii sufficiently large to encompass the cosmic baryon fraction, but it is unclear where or even if this occurs in galaxies. Thus, while the $\Lambda$CDM prediction for the full RAR can be tested only by applying potentially restrictive priors on galaxy formation effects (Di Cintio & Lelli 2016; Desmond 2017; Paranjape & Sheth 2021), a more direct route towards informing the dark matter-modified gravity debate is to test the MOND prediction, specifically the limiting behaviour at $g \ll g_0$ and $g \gg g_0$, the small intrinsic scatter and the lack of residual correlations.

Here, we focus on the asymptotic behaviour. This can be assessed to some extent by fitting a functional form with free power-law slopes at both ends (Lelli et al. 2017), but this assumes that the slope tends to a constant at each end and restricts to a specific part of the functional parameter space for which this is the case. These are in question when assessing the accuracy of the MOND prescription. A fully satisfactory fit should therefore make no such assumptions, eliminating potential confirmation bias and testing without any priors the assertion that the RAR implies no dynamically relevant dark matter at high $g$ and the deep-MOND limit at low $g$. We accomplish this here by means of a novel regression algorithm dubbed *Exhaustive Symbolic Regression* (ESR; Bartlett, Desmond & Ferreira 2022b), and hence assess the degree to which the RAR supports the tenets of MOND. Within the MOND paradigm, this method also enables optimization of the ‘interpolating function’ (IF) $g_{\text{obs}} = F(g_{\text{bar}})$ between the two stipulated limits.

The structure of the paper is as follows. In Section 2, we describe the RAR data that we use, and in Section 3, our algorithm for generating functions and assessing their aptitude for describing the data. Section 4 presents the results. In Section 5, we discuss the broader ramifications, potential remaining uncertainties and ways in which the programme could be furthered in the future. Section 6 concludes; full details on ESR are given in the companion paper Bartlett et al. (2022b). Units not explicitly given are $10^{-10} \text{ m s}^{-2}$, and all logarithms are natural.

## 2 OBSERVATIONAL DATA

We use the SPARC data set (Lelli, McGaugh & Schombert 2016),<fnref n="1" /> a compilation of 175 rotation curves from the literature combined with Spitzer 3.6 µm photometry. We apply the same quality cuts as the RAR study of Lelli et al. (2017), removing galaxies with quality flag 3 (indicating large asymmetries, non-circular motions and/or offsets between stellar and HI distributions), and those with inclinations $i < 30^\circ$ 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.

## 3 METHOD

We describe our method for generating and assessing trial functions in Section 3.1, and our likelihood function in Section 3.2. In Section 3.3, we outline our criteria for assessing whether a function displays MOND-like behaviour.

### 3.1 Exhaustive symbolic regression

While algorithms for symbolic regression (SR) – the search for good functional descriptions of a data set – are becoming mature, they

remain fallible (La Cava et al. 2021). Unless the generating function of the data is known at the outset (in which case SR is not required), it is not possible to determine whether any SR algorithm has uncovered the best function. This motivated us to develop ESR which, given a set of basis functions, produces and evaluates every possible function up to a given complexity of equation, defined here as the number of nodes in its tree representation. This enables a brute-force solution to relatively simple problems and provides a touchstone for assessing the results of stochastic algorithms at higher complexity. As shown in detail in Bartlett et al. (2022b), stochastic searches regularly fail to the best functions at even moderate complexity ~6, so that we would not be confidant of obtaining the functional form of the RAR through any algorithm besides ESR.

Presented in full in the companion paper, ESR has two main steps: i) generating, and optimising the parameters of, all functions up to a given complexity, and ii) ranking these functions using an information-theoretic metric combining accuracy and simplicity. For part i, the steps are:

(i) Generate all possible trees containing a given number of nodes (equal to the complexity of functions considered).

(ii) Generate the complete set of such functions by decorating these trees with all permutations of the operators from the operator list specified in advance, utilising the constraints on the arity of the operator that can occupy a given node.

(iii) Simplify the functions and remove duplicates. Variants of the same function (e.g. $x(x + \theta_0)$ and $x^2 + \theta_0 x$) are however retained as these may have different model complexities (used in step ii, below). For each unique function the variant is retained that minimizes this.

(iv) Determine the values of the free parameters appearing in the functions that maximize the likelihood of the data (see Section 3.2).

(v) Repeat for all complexities under consideration.

The only degrees of freedom in this procedure are the maximum complexity considered (here set at 9 as higher complexity is computationally prohibitive) and the set of operators of which the functions are composed. Here we choose:<fnref n="2" />

(i) **Nullary**: $g_{\text{bar}}, \theta$

(ii) **Unary**: exp, sqrt, square, inv

(iii) **Binary**: $+, -, *, /, pow$

where $\theta$ is a free parameter. We implicitly take the absolute value of the argument of any square root or power.

The result of this procedure is a list of all functions up to the maximum complexity (of which there are $2.24 \times 10^7$), along with the parameter values that maximize the likelihood of the RAR data. As in regular regression, using the maximum likelihood as the model selection criterion would favour overfitting, whereby a function fits the data near-perfectly but generalizes or extrapolates poorly. To remedy this, SR typically uses two-objective optimization, where the second objective is the ‘simplicity’ of the function. In the absence of a metric for trading accuracy (the first objective) with complexity, optimal functions form a ‘Pareto front’ where accuracy cannot be increased without reducing simplicity and vice versa. Simplicity has been defined analogously to model complexity (the number of nodes in the tree representation; e.g. in PySR; Cranmer et al. 2020), among others, but such definitions are typically

<footnote n="1">¹http://astroweb.cwru.edu/SPARC/</footnote>

<footnote n="2">²Many of these operators can manifestly be constructed from combinations of others. We include these repetitions to simplify the function trees, allowing a greater range of expressions up to the maximum complexity.</footnote>
arbitrary and thus compromise the objectivity of the regression results.

To remedy this, part (ii) of ESR implements the minimum description length principle (MDL; Rissanen 1978; Grunwald 2007; Grünwald & Roos 2019) as a model selection criterion, which has an information-theoretic motivation and provides a natural framework for making commensurable the two objectives. MDL states that functions are preferred to the extent that they compress the data, i.e. minimize the number of bits required to communicate the data with the aid of the function. We implement this with a two-step code in which the description length (also called codelength) is comprised of a component describing the function and a component describing the residuals of the data around the function’s expectation. We use the Shannon–Fano coding scheme for the latter (Cover & Thomas 1991), and for the former include contributions both from the structure of the function (penalising those employing more operators) and from the free parameters (penalising more parameters, especially ones that must be specified to high precision to achieve a high likelihood). The overall codelength of the compressed data, $L(D)$, is derived in section 3 of Bartlett et al. (2022b):

$$ 
\begin{aligned}
L(D) &= L(D|H) + L(H) \\
&= -\log(\hat{\mathcal{L}}) + k \log(n) - \frac{p}{2} \log(3) \\
&\quad + \sum_{i}^{p} \left( \frac{1}{2} \log(\hat{I}_{ii}) + \log(|\hat{\theta}_{i}|) \right) + \sum_{j} \log(c_j),
\end{aligned} 
\quad (1) 
$$

where $L$ is the description length, $D$ the data set, $H$ the hypothesis (i.e. function in question), $\mathcal{L}$ the likelihood, $\theta$ a free parameter of the function, $k$ the number of nodes in the function’s tree representation, $n$ the number of unique operators involved, $p$ the total number of free parameters, $I$ the Fisher information matrix of the parameters and $c_j$ any constant natural numbers generated by simplifications. A hat denotes evaluation at the maximum-likelihood point. With all logarithms natural, this is the number of nats required to communicate the data with the aid of the function. $L(D)$ supports a probabilistic interpretation over function space that generalizes the likelihood: the relative probability of a function is $\exp(-L(D))$ (Grunwald 2007).

The structure of the function alone determines the klog(n) term, but the remaining terms require the free parameters to be numerically optimized to maximize the likelihood (which we use interchangeably with minimising the loss).<fnref n="3" /> We now describe our choice of likelihood for the RAR data.

## 3.2 Loss function

As is typical (e.g. Lelli et al. 2017), we assume that $g_{\text{bar}}$, $g_{\text{obs}}$ and their uncertainties are uncorrelated across the data set. We further assume that the true $g_{\text{bar}}$ and $g_{\text{obs}}$ values, denoted $g_{\text{bar}}^{\text{t}}$ and $g_{\text{obs}}^{\text{t}}$, generate the observed values with lognormal probability distributions centred at the true values with widths given by their uncertainties $\delta g_{\text{bar}}$ and $\delta g_{\text{obs}}$. Following Lelli et al. (2017), we fix the mass-to-light ratios $\Upsilon_{\text{gas}} = 1.33$, $\Upsilon_{\text{disc}} = 0.5$, and $\Upsilon_{\text{bulge}} = 0.7$, and assign them 10, 25, and 25 per cent uncertainties, respectively, summing these in quadrature to estimate $\delta V_{\text{bar}}$ and hence $\delta g_{\text{bar}}$ (assuming no uncertainty in radial position). We likewise assume the uncertainties on distance $D$ and

inclination $i$ to be statistical and hence sum their contributions in quadrature with the quoted statistical uncertainty on $V_{\text{obs}}$ according to Lelli et al. (2017, equation 2) to estimate $\delta g_{\text{bar}}$.

The likelihood of an observation given the function in question, $f(g_{\text{bar}}^{\text{t}})$, is then:

$$ 
\begin{aligned}
\mathcal{L}(\log(g_{\text{obs}})) &= \int_{-\infty}^{\infty} \mathcal{L}(\log(g_{\text{obs}})|\log(g_{\text{bar}}^{\text{t}})) \mathcal{L}(\log(g_{\text{bar}}^{\text{t}})) \mathrm{d} \log(g_{\text{bar}}^{\text{t}}) \\
&= \frac{1}{2\pi \delta \log(g_{\text{bar}}) \delta \log(g_{\text{obs}})} \\
&\quad \times \int_{-\infty}^{\infty} \exp \left( - \frac{(\log(g_{\text{obs}}) - \log(f(g_{\text{bar}}^{\text{t}})))^2}{2 \delta \log(g_{\text{obs}})^2} \right) \\
&\quad \times \exp \left( - \frac{(\log(g_{\text{bar}}^{\text{t}}) - \log(g_{\text{bar}}))^2}{2 \delta \log(g_{\text{bar}})^2} \right) \mathrm{d} \log(g_{\text{bar}}^{\text{t}}) \\
&\approx \frac{1}{\sqrt{2\pi \sigma_{\text{tot}}^2}} \exp \left( - \frac{(\log(g_{\text{obs}}) - \log(f(g_{\text{bar}}^{\text{t}})))^2}{2\sigma_{\text{tot}}^2} \right),
\end{aligned} 
\quad (2) 
$$

where

$$ \sigma_{\text{tot}}^2 \equiv \delta \log(g_{\text{obs}})^2 + \left( \frac{\mathrm{d} \log(f(g_{\text{bar}}))}{\mathrm{d} \log(g_{\text{bar}})} \right)^2 \delta \log(g_{\text{bar}})^2. \quad (3) $$

This is derived by keeping the leading-order term in the Taylor expansion of $\log(f(g_{\text{bar}}^{\text{t}}))$ around $\log(f(g_{\text{bar}}))$, and therefore assumes this to be small relative to the rate of change of $f$. An advantage of working in $\log(g_{\text{bar}}) - \log(g_{\text{obs}})$ space as opposed to $g_{\text{bar}} - g_{\text{obs}}$ is that, the RAR being roughly a product of power-laws, this minimizes the error due to the first-order approximation. We find this to be good for all of the best functions. The likelihood is then the product over all data points. We discuss in Section A, the limitations of this likelihood model and how it could be improved.

## 3.3 Assessing MOND

The core of the MOND paradigm is that $g_{\text{obs}} = g_{\text{bar}}$ at $g \gg g_0$ and $g_{\text{obs}} = \sqrt{g_{\text{bar}}g_0}$ at $g \ll g_0$ (Milgrom 1983a, b, c). This implies

$$ s \equiv \frac{\mathrm{d} \log(g_{\text{obs}})}{\mathrm{d} \log(g_{\text{bar}})} = \begin{cases} 1, & g_{\text{bar}} \to \infty \\ 1/2, & g_{\text{bar}} \to 0. \end{cases} \quad (4) $$

Common choices for the IF covering the intermediate region $g \approx g_0$ include

(i) **'Simple':** $g_{\text{obs}} = g_{\text{bar}}/2 + \sqrt{g_{\text{bar}}^2/4 + g_{\text{bar}}g_0}$ (Famaey & Binney 2005),

(ii) **'Standard':** $g_{\text{obs}} = \frac{1}{\sqrt{2}} \sqrt{g_{\text{bar}}^2 + \sqrt{g_{\text{bar}}^2(g_{\text{bar}}^2 + 4g_0^2)}}$ (Milgrom 1983c),

(iii) **'RAR':** $g_{\text{obs}} = g_{\text{bar}}/(1 - \exp(-\sqrt{g_{\text{bar}}/g_0}))$ (Lelli et al. 2017).

Fig. 1 plots these functions on top of the RAR data for the best-fitting values on the SPARC data (shown in the lower rows of Table 1 in Section 4.1). The Simple and RAR IFs are distinguished from the Standard IF principally by a more gradual transition between the Newtonian and deep-MOND regimes, although the Standard IF also prefers a significantly higher value of $g_0$. While the basic MOND framework is not committed to any particular IF, it is committed to equation (4) providing an optimal description of the data. Our assessment of the theory will therefore be based on the extent to which the best functions (those with lowest description length) conform to these limits: any function that does so, in addition to possessing a coefficient of proportionality of unity in $g_{\text{obs}} \propto g_{\text{bar}}$ at high $g_{\text{bar}}$, may be considered a new MOND IF (The low-$g_{\text{bar}}$

<footnote n="3">³The $p \log(3)/2$ term is only affected by the numerical optimization if any parameters are set to 0 due to their maximum-likelihood values being less than 1 precision unit from 0 (see section 3 of Bartlett et al. 2022b).</footnote>
arbitrary and thus compromise the objectivity of the regression results.

To remedy this, part (ii) of ESR implements the minimum description length principle (MDL; Rissanen 1978; Grunwald 2007; Grünwald & Roos 2019) as a model selection criterion, which has an information-theoretic motivation and provides a natural framework for making commensurable the two objectives. MDL states that functions are preferred to the extent that they compress the data, i.e. minimize the number of bits required to communicate the data with the aid of the function. We implement this with a two-step code in which the description length (also called codelength) is comprised of a component describing the function and a component describing the residuals of the data around the function’s expectation. We use the Shannon–Fano coding scheme for the latter (Cover & Thomas 1991), and for the former include contributions both from the structure of the function (penalising those employing more operators) and from the free parameters (penalising more parameters, especially ones that must be specified to high precision to achieve a high likelihood). The overall codelength of the compressed data, $L(D)$, is derived in section 3 of Bartlett et al. (2022b):

$$ 
\begin{aligned}
L(D) &= L(D|H) + L(H) \\
&= -\log(\hat{\mathcal{L}}) + k \log(n) - \frac{p}{2} \log(3) \\
&\quad + \sum_{i}^{p} \left( \frac{1}{2} \log(\hat{I}_{ii}) + \log(|\hat{\theta}_{i}|) \right) + \sum_{j} \log(c_j),
\end{aligned} 
\quad (1) 
$$

where $L$ is the description length, $D$ the data set, $H$ the hypothesis (i.e. function in question), $\mathcal{L}$ the likelihood, $\theta$ a free parameter of the function, $k$ the number of nodes in the function’s tree representation, $n$ the number of unique operators involved, $p$ the total number of free parameters, $I$ the Fisher information matrix of the parameters and $c_j$ any constant natural numbers generated by simplifications. A hat denotes evaluation at the maximum-likelihood point. With all logarithms natural, this is the number of nats required to communicate the data with the aid of the function. $L(D)$ supports a probabilistic interpretation over function space that generalizes the likelihood: the relative probability of a function is $\exp(-L(D))$ (Grunwald 2007).

The structure of the function alone determines the $k\log(n)$ term, but the remaining terms require the free parameters to be numerically optimized to maximize the likelihood (which we use interchangeably with minimising the loss).<fnref n="3" /> We now describe our choice of likelihood for the RAR data.

## 3.2 Loss function

As is typical (e.g. Lelli et al. 2017), we assume that $g_{\text{bar}}$, $g_{\text{obs}}$ and their uncertainties are uncorrelated across the data set. We further assume that the true $g_{\text{bar}}$ and $g_{\text{obs}}$ values, denoted $g_{\text{bar}}^{\text{t}}$ and $g_{\text{obs}}^{\text{t}}$, generate the observed values with lognormal probability distributions centred at the true values with widths given by their uncertainties $\delta g_{\text{bar}}$ and $\delta g_{\text{obs}}$. Following Lelli et al. (2017), we fix the mass-to-light ratios $\Upsilon_{\text{gas}} = 1.33$, $\Upsilon_{\text{disc}} = 0.5$, and $\Upsilon_{\text{bulge}} = 0.7$, and assign them 10, 25, and 25 per cent uncertainties, respectively, summing these in quadrature to estimate $\delta V_{\text{bar}}$ and hence $\delta g_{\text{bar}}$ (assuming no uncertainty in radial position). We likewise assume the uncertainties on distance $D$ and

inclination $i$ to be statistical and hence sum their contributions in quadrature with the quoted statistical uncertainty on $V_{\text{obs}}$ according to Lelli et al. (2017, equation 2) to estimate $\delta g_{\text{bar}}$.

The likelihood of an observation given the function in question, $f(g_{\text{bar}}^{\text{t}})$, is then:

$$ 
\begin{aligned}
\mathcal{L}(\log(g_{\text{obs}})) &= \int_{-\infty}^{\infty} \mathcal{L}(\log(g_{\text{obs}})|\log(g_{\text{bar}}^{\text{t}})) \mathcal{L}(\log(g_{\text{bar}}^{\text{t}})) \mathrm{d}\log(g_{\text{bar}}^{\text{t}}) \\
&= \frac{1}{2\pi \delta \log(g_{\text{bar}}) \delta \log(g_{\text{obs}})} \\
&\quad \times \int_{-\infty}^{\infty} \exp \left( -\frac{(\log(g_{\text{obs}}) - \log(f(g_{\text{bar}}^{\text{t}})))^2}{2 \delta \log(g_{\text{obs}})^2} \right) \\
&\quad \times \exp \left( -\frac{(\log(g_{\text{bar}}^{\text{t}}) - \log(g_{\text{bar}}))^2}{2 \delta \log(g_{\text{bar}})^2} \right) \mathrm{d}\log(g_{\text{bar}}^{\text{t}}) \\
&\approx \frac{1}{\sqrt{2\pi\sigma_{\text{tot}}^2}} \exp \left( -\frac{(\log(g_{\text{obs}}) - \log(f(g_{\text{bar}}^{\text{t}})))^2}{2\sigma_{\text{tot}}^2} \right),
\end{aligned} 
\quad (2) 
$$

where

$$ \sigma_{\text{tot}}^2 \equiv \delta \log(g_{\text{obs}})^2 + \left( \frac{\mathrm{d} \log(f(g_{\text{bar}}))}{\mathrm{d} \log(g_{\text{bar}})} \right)^2 \delta \log(g_{\text{bar}})^2. \quad (3) $$

This is derived by keeping the leading-order term in the Taylor expansion of $\log(f(g_{\text{bar}}^{\text{t}}))$ around $\log(f(g_{\text{bar}}))$, and therefore assumes this to be small relative to the rate of change of $f$. An advantage of working in $\log(g_{\text{bar}}) - \log(g_{\text{obs}})$ space as opposed to $g_{\text{bar}} - g_{\text{obs}}$ is that, the RAR being roughly a product of power-laws, this minimizes the error due to the first-order approximation. We find this to be good for all of the best functions. The likelihood is then the product over all data points. We discuss in Section A, the limitations of this likelihood model and how it could be improved.

## 3.3 Assessing MOND

The core of the MOND paradigm is that $g_{\text{obs}} = g_{\text{bar}}$ at $g \gg g_0$ and $g_{\text{obs}} = \sqrt{g_{\text{bar}}g_0}$ at $g \ll g_0$ (Milgrom 1983a, b, c). This implies

$$ s \equiv \frac{\mathrm{d} \log(g_{\text{obs}})}{\mathrm{d} \log(g_{\text{bar}})} = \begin{cases} 1, & g_{\text{bar}} \to \infty \\ 1/2, & g_{\text{bar}} \to 0. \end{cases} \quad (4) $$

Common choices for the IF covering the intermediate region $g \approx g_0$ include

(i) **'Simple':** $g_{\text{obs}} = g_{\text{bar}}/2 + \sqrt{g_{\text{bar}}^2/4 + g_{\text{bar}}g_0}$ (Famaey & Binney 2005),

(ii) **'Standard':** $g_{\text{obs}} = \frac{1}{\sqrt{2}} \sqrt{g_{\text{bar}}^2 + \sqrt{g_{\text{bar}}^2(g_{\text{bar}}^2 + 4g_0^2)}}$ (Milgrom 1983c),

(iii) **'RAR':** $g_{\text{obs}} = g_{\text{bar}}/(1 - \exp(-\sqrt{g_{\text{bar}}/g_0}))$ (Lelli et al. 2017).

Fig. 1 plots these functions on top of the RAR data for the best-fitting values on the SPARC data (shown in the lower rows of Table 1 in Section 4.1). The Simple and RAR IFs are distinguished from the Standard IF principally by a more gradual transition between the Newtonian and deep-MOND regimes, although the Standard IF also prefers a significantly higher value of $g_0$. While the basic MOND framework is not committed to any particular IF, it is committed to equation (4) providing an optimal description of the data. Our assessment of the theory will therefore be based on the extent to which the best functions (those with lowest description length) conform to these limits: any function that does so, in addition to possessing a coefficient of proportionality of unity in $g_{\text{obs}} \propto g_{\text{bar}}$ at high $g_{\text{bar}}$, may be considered a new MOND IF (The low-$g_{\text{bar}}$

<footnote n="3">³The $p \log(3)/2$ term is only affected by the numerical optimization if any parameters are set to 0 due to their maximum-likelihood values being less than 1 precision unit from 0 (see section 3 of Bartlett et al. 2022b).</footnote>
coefficient of proportionality is $\sqrt{g_0}$ but $g_0$ is unknown a priori, so this does not supply an additional requirement.). Following Lelli et al. (2017), we will also consider a double power law fit:

$$g_{\text{obs}} = \theta_1 \left(1 + \frac{g_{\text{bar}}}{\theta_0}\right)^{\theta_2 - \theta_3} \left(\frac{g_{\text{bar}}}{\theta_0}\right)^{\theta_3}, \quad (5)$$

which has limiting logarithmic slopes of $\theta_3$ and $\theta_2$, and plot the best-fit in Fig. 1. We define $s_- \equiv \lim_{g_{\text{bar}} \to 0} s$ and $s_+ \equiv \lim_{g_{\text{bar}} \to \infty} s$.

The MOND interpretation of the RAR is complicated by the possibility of the external field effect (EFE), a breakdown of the strong equivalence principle due to the non-linear, acceleration-based modification to Newtonian mechanics (Milgrom 1983a). The EFE implies that otherwise identical galaxies in different external gravitational fields have different dynamics, which is a function of the external field strength $g_{\text{ex}}$ relative to $g_0$ and the internal field $g_{\text{in}}$. In the quasi-Newtonian regime $g_{\text{in}} < g_{\text{ex}} < g_0$, Kepler's laws are recovered with dynamical masses scaled by $g_{\text{ex}}/g_0$, while in the external field-dominated regime $g_{\text{in}} < g_0 < g_{\text{ex}}$, Newtonian mechanics are fully recovered (Famaey & McGaugh 2012). This steepens the RAR at low $g_{\text{bar}}$.

The precise effect of the EFE is difficult to calculate in general because it depends both on the underlying MOND theory, and on a galaxy's morphology and orientation with respect to the external field direction. The most sophisticated fitting functions to MOND simulations are currently to be found in Zonoozi et al. (2021) for QUMOND (Milgrom 2010) and Chae & Milgrom (2022) for AQUAL (Bekenstein & Milgrom 1984). Chae et al. (2022) tested these expectations by fitting the SPARC rotation curves for the average external field strength, finding this to be in good agreement with independent estimates based on the baryonic mass surrounding the SPARC galaxies (Chae et al. 2021) for AQUAL, but less so for QUMOND. We will therefore use AQUAL to explore the effect of the EFE on the expected low-$g_{\text{bar}}$ slope and functional form more generally. Chae & Milgrom (2022) eq. 15 gives

$$g_{\text{obs}} = g_{\text{bar}} \left( \frac{1}{2} + \left( \frac{1}{4} + \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{\frac{1}{2}} \right) \times \left( 1 + \tanh \left( \frac{1.1e_N}{g_{\text{bar}}/g_0} \right)^{1.2} \times \left( -\frac{1}{3} \right) \times \frac{\left( \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right) \left( \frac{1}{4} + \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{-\frac{1}{2}}}{1 + \left( \frac{1}{2} + 2 \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{\frac{1}{2}}} \right). \quad (6)$$

This allows for variable disc thickness and scale length, and is azimuthally averaged to reduce sensitivity to the orientation of the field relative to the disc axis. It recovers the Simple IF as $e_N \equiv g_{\text{ex}}/g_0 \to 0$ and hence we refer to it as 'Simple IF + EFE'.

Note that, while some form of the EFE is generically predicted by MOND in modified inertia formulations, it may be very different (e.g. a function of the entire past trajectory of an object) or negligible (Milgrom 2011). While there is evidence for the EFE in many systems (McGaugh & Milgrom 2013; Haghi et al. 2019; Chae et al. 2020b), in others it appears conspicuously absent (Hernandez et al. 2019; Freundlich et al. 2022). The black curve in Fig. 1 shows the best fit to the data using equation (6).

<figure><img src="image_1.png" /><figcaption>Figure 1. The Simple, Standard and RAR IFs, a double power law, and Simple IF with a global external field strength in AQUAL, overlaid on the SPARC data (blue points). The parameters are set to their maximum-likelihood values shown in Table 1. The dashed black line shows the one-to-one relation (Newtonian limit) and the cross in the lower right shows the average uncertainty size.</figcaption></figure>

## 3.4 Mock data generation

To shed light on the significance of our results, we apply ESR also to two stacks of mock data sets. We generate each mock data set using exactly the same number of points as the SPARC data, and with identical $\log g_{\text{bar}}$, $\delta \log g_{\text{bar}}$, and $\delta \log g_{\text{obs}}$ values, but with $\log g_{\text{obs}}$ generated using a MONDian function. This assumes $g_{\text{bar}}^t$ equal to the SPARC $g_{\text{bar}}$ (the maximum a priori estimate), $\log g_{\text{bar}}$ for each mock realization drawn from $N(\log g_{\text{bar}}^t, \delta \log g_{\text{bar}})$, and $\log g_{\text{obs}}$ from $N(F(\log g_{\text{bar}}^t), \delta \log g_{\text{obs}})$. To reduce the impact of noise in the mock data, we apply ESR to a stack of 10 independent realizations.<fnref n="4" /> The only terms in the description length that depend on the data set size are $\log(\mathcal{L})$ and the Fisher matrix $I$, both of which scale linearly. To make the results compatible with the real data we therefore divide these terms by 10.

The two mock data set stacks differ in the generating function $\mathcal{F}$. For the first, we use the RAR IF with the best-fitting value on the data $g_0 = 1.127$ (see Section 4). This function is already known to describe the RAR well (Lelli et al. 2017), and has low enough complexity to be included (as a special case of a more general function, see below) in our function list. Since equation (4) is satisfied by construction in this case, evaluating it on the best functions from ESR will address the question of whether the dynamic range of the data is sufficiently high – and the uncertainties sufficiently low – to pick out unambiguously a correctly MONDian solution, as only in this case could one expect to obtain such behaviour for the real data were it generated by MOND.

The second stack is created using equation (6). We adopt $g_0 = 1.2$ and $\langle g_{\text{ex}} \rangle = 1.2 \times 10^{-2}$ ($e_N = 0.01$), corresponding roughly to maximal clustering of unobserved baryons (as expected in a

<footnote n="4">⁴The reason not to do more is that the time required for parameter optimization scales with the number of data points, and this is already expensive at complexity 9. We assess convergence by fitting the complexity 4–6 equations to an independent stack of 10 realizations. We find that the order of equations when sorted by $L(D)$ is identical between the two stacks, and that the difference between the description lengths of particular equations falls with complexity, from ~40 at complexity 4 to ~5 at complexity 6. This implies that the best equations, mostly at complexity 8 and 9, are not sensitive to the random number generation.</footnote>
coefficient of proportionality is $\sqrt{g_0}$ but $g_0$ is unknown a priori, so this does not supply an additional requirement.). Following Lelli et al. (2017), we will also consider a double power law fit:

$$g_{\text{obs}} = \theta_1 \left(1 + \frac{g_{\text{bar}}}{\theta_0}\right)^{\theta_2 - \theta_3} \left(\frac{g_{\text{bar}}}{\theta_0}\right)^{\theta_3}, \quad (5)$$

which has limiting logarithmic slopes of $\theta_3$ and $\theta_2$, and plot the best-fit in Fig. 1. We define $s_- \equiv \lim_{g_{\text{bar}} \to 0} s$ and $s_+ \equiv \lim_{g_{\text{bar}} \to \infty} s$.

The MOND interpretation of the RAR is complicated by the possibility of the external field effect (EFE), a breakdown of the strong equivalence principle due to the non-linear, acceleration-based modification to Newtonian mechanics (Milgrom 1983a). The EFE implies that otherwise identical galaxies in different external gravitational fields have different dynamics, which is a function of the external field strength $g_{\text{ex}}$ relative to $g_0$ and the internal field $g_{\text{in}}$. In the quasi-Newtonian regime $g_{\text{in}} < g_{\text{ex}} < g_0$, Kepler's laws are recovered with dynamical masses scaled by $g_{\text{ex}}/g_0$, while in the external field-dominated regime $g_{\text{in}} < g_0 < g_{\text{ex}}$, Newtonian mechanics are fully recovered (Famaey & McGaugh 2012). This steepens the RAR at low $g_{\text{bar}}$.

The precise effect of the EFE is difficult to calculate in general because it depends both on the underlying MOND theory, and on a galaxy's morphology and orientation with respect to the external field direction. The most sophisticated fitting functions to MOND simulations are currently to be found in Zonoozi et al. (2021) for QUMOND (Milgrom 2010) and Chae & Milgrom (2022) for AQUAL (Bekenstein & Milgrom 1984). Chae et al. (2022) tested these expectations by fitting the SPARC rotation curves for the average external field strength, finding this to be in good agreement with independent estimates based on the baryonic mass surrounding the SPARC galaxies (Chae et al. 2021) for AQUAL, but less so for QUMOND. We will therefore use AQUAL to explore the effect of the EFE on the expected low-$g_{\text{bar}}$ slope and functional form more generally. Chae & Milgrom (2022) eq. 15 gives

$$g_{\text{obs}} = g_{\text{bar}} \left( \frac{1}{2} + \left( \frac{1}{4} + \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{\frac{1}{2}} \right) \times \left( 1 + \tanh \left( \frac{1.1e_N}{g_{\text{bar}}/g_0} \right)^{1.2} \times \left( -\frac{1}{3} \right) \times \frac{\left( \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right) \left( \frac{1}{4} + \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{-\frac{1}{2}}}{1 + \left( \frac{1}{2} + 2 \left( \left( \frac{g_{\text{bar}}}{g_0} \right)^2 + (1.1e_N)^2 \right)^{-\frac{1}{2}} \right)^{\frac{1}{2}}} \right). \quad (6)$$

This allows for variable disc thickness and scale length, and is azimuthally averaged to reduce sensitivity to the orientation of the field relative to the disc axis. It recovers the Simple IF as $e_N \equiv g_{\text{ex}}/g_0 \to 0$ and hence we refer to it as 'Simple IF + EFE'.

Note that, while some form of the EFE is generically predicted by MOND in modified inertia formulations, it may be very different (e.g. a function of the entire past trajectory of an object) or negligible (Milgrom 2011). While there is evidence for the EFE in many systems (McGaugh & Milgrom 2013; Haghi et al. 2019; Chae et al. 2020b), in others it appears conspicuously absent (Hernandez et al. 2019; Freundlich et al. 2022). The black curve in Fig. 1 shows the best fit to the data using equation (6).

<figure><img src="image_2.png" /><figcaption>Figure 1. The Simple, Standard and RAR IFs, a double power law, and Simple IF with a global external field strength in AQUAL, overlaid on the SPARC data (blue points). The parameters are set to their maximum-likelihood values shown in Table 1. The dashed black line shows the one-to-one relation (Newtonian limit) and the cross in the lower right shows the average uncertainty size.</figcaption></figure>

## 3.4 Mock data generation

To shed light on the significance of our results, we apply ESR also to two stacks of mock data sets. We generate each mock data set using exactly the same number of points as the SPARC data, and with identical $\log g_{\text{bar}}$, $\delta \log g_{\text{bar}}$, and $\delta \log g_{\text{obs}}$ values, but with $\log g_{\text{obs}}$ generated using a MONDian function. This assumes $g_{\text{bar}}^t$ equal to the SPARC $g_{\text{bar}}$ (the maximum a priori estimate), $\log g_{\text{bar}}$ for each mock realization drawn from $N(\log g_{\text{bar}}^t, \delta \log g_{\text{bar}})$, and $\log g_{\text{obs}}$ from $N(F(\log g_{\text{bar}}^t), \delta \log g_{\text{obs}})$. To reduce the impact of noise in the mock data, we apply ESR to a stack of 10 independent realizations.<fnref n="4" /> The only terms in the description length that depend on the data set size are $\log(L)$ and the Fisher matrix $I$, both of which scale linearly. To make the results compatible with the real data we therefore divide these terms by 10.

The two mock data set stacks differ in the generating function $F$. For the first, we use the RAR IF with the best-fitting value on the data $g_0 = 1.127$ (see Section 4). This function is already known to describe the RAR well (Lelli et al. 2017), and has low enough complexity to be included (as a special case of a more general function, see below) in our function list. Since equation (4) is satisfied by construction in this case, evaluating it on the best functions from ESR will address the question of whether the dynamic range of the data is sufficiently high – and the uncertainties sufficiently low – to pick out unambiguously a correctly MONDian solution, as only in this case could one expect to obtain such behaviour for the real data were it generated by MOND.

The second stack is created using equation (6). We adopt $g_0 = 1.2$ and $\langle g_{\text{ex}} \rangle = 1.2 \times 10^{-2}$ ($e_N = 0.01$), corresponding roughly to maximal clustering of unobserved baryons (as expected in a

<footnote n="4">⁴The reason not to do more is that the time required for parameter optimization scales with the number of data points, and this is already expensive at complexity 9. We assess convergence by fitting the complexity 4–6 equations to an independent stack of 10 realizations. We find that the order of equations when sorted by $L(D)$ is identical between the two stacks, and that the difference between the description lengths of particular equations falls with complexity, from ~40 at complexity 4 to ~5 at complexity 6. This implies that the best equations, mostly at complexity 8 and 9, are not sensitive to the random number generation.</footnote>
<table><caption>**Table 1.** Top functions found by ESR applied to the SPARC data, ranked by total description length, compared to four MOND IFs, and a double power law below the horizontal line. $x \equiv g_{bar}/10^{-10} m s^{-2}$. We include also the 'generalized RAR IF' at rank 17, although this does not have a deep-MOND limit. The IFs and double power law are not produced explicitly by our implementation of ESR and hence their ranks are unknown, although they are clearly worse than the best ESR functions in both description length and likelihood. For the Simple, RAR and Standard IFs the parameter is $g_0$, while for the 'Simple IF + EFE' (equation 6), the first is $g_0$ and the second $e_N$.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th rowspan="2">Comp.</th><th rowspan="2">P(f)</th><th colspan="4">Parameters</th><th rowspan="2">Resid.<sup>1</sup></th><th colspan="2">Description length</th><th rowspan="2">Total</th></tr><tr><th>θ<sub>0</sub></th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th><th>θ<sub>3</sub></th><th>Func.<sup>2</sup></th><th>Param.<sup>3</sup></th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> (|θ<sub>1</sub> + x|<sup>2</sup> + x)</td><td>9</td><td>9.3 × 10<sup>-1</sup></td><td>0.84</td><td>-0.02</td><td>0.38</td><td>-</td><td>-1279.1</td><td>14.5</td><td>14.0</td><td>-1250.6</td></tr><tr><td>2</td><td>|θ<sub>1</sub>|<sup>x</sup> + θ<sub>0</sub>|<sup>2</sup> + x</td><td>9</td><td>6.4 × 10<sup>-2</sup></td><td>-0.99</td><td>0.64</td><td>0.36</td><td>-</td><td>-1279.9</td><td>12.5</td><td>19.6</td><td>-1247.9</td></tr><tr><td>3</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub> - x|<sup>2</sup> - θ<sub>3</sub></sup></td><td>9</td><td>2.0 × 10<sup>-3</sup></td><td>-1.4 × 10<sup>2</sup></td><td>0.02</td><td>0.14</td><td>0.89</td><td>-1276.4</td><td>12.5</td><td>19.5</td><td>-1244.4</td></tr><tr><td>4</td><td>|θ<sub>0</sub>(θ<sub>1</sub> + x)|<sup>2</sup> + x</td><td>9</td><td>1.4 × 10<sup>-4</sup></td><td>0.35</td><td>-0.02</td><td>0.34</td><td>-</td><td>-1268.9</td><td>14.5</td><td>12.7</td><td>-1241.7</td></tr><tr><td>5</td><td>|θ<sub>0</sub> - |θ<sub>1</sub> - x|<sup>2</sup>|<sup>3</sup></td><td>9</td><td>1.0 × 10<sup>-5</sup></td><td>-0.30</td><td>0.02</td><td>0.42</td><td>2.14</td><td>-1271.1</td><td>12.5</td><td>19.5</td><td>-1239.1</td></tr><tr><td>6</td><td>√x exp (|θ<sub>0</sub>+x|<sup>1</sup>/2)</td><td>9</td><td>1.5 × 10<sup>-9</sup></td><td>-0.02</td><td>0.36</td><td>-</td><td>-</td><td>-1257.9</td><td>17.5</td><td>10.0</td><td>-1230.3</td></tr><tr><td>7</td><td>(|θ<sub>0</sub>|<sup>x</sup>/x) + x</td><td>9</td><td>2.4 × 10<sup>-10</sup></td><td>1.87</td><td>-0.52</td><td>-</td><td>-</td><td>-1250.6</td><td>14.5</td><td>7.6</td><td>-1228.5</td></tr><tr><td>8</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>8</td><td>1.8 × 10<sup>-10</sup></td><td>-1.8 × 10<sup>-3</sup></td><td>0.72</td><td>-</td><td>-</td><td>-1245.6</td><td>12.9</td><td>4.5</td><td>-1228.2</td></tr><tr><td>9</td><td>|θ<sub>0</sub> + 1/√x|<sup>θ<sub>1</sub></sup></td><td>8</td><td>9.6 × 10<sup>-11</sup></td><td>-0.22</td><td>-2.14</td><td>-</td><td>-</td><td>-1251.1</td><td>14.3</td><td>9.2</td><td>-1227.6</td></tr><tr><td>10</td><td>(√x + 1/√x)<sup>θ<sub>0</sub></sup> + x</td><td>9</td><td>8.2 × 10<sup>-11</sup></td><td>-0.53</td><td>-</td><td>-</td><td>-</td><td>-1248.3</td><td>16.1</td><td>4.8</td><td>-1227.4</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>17</td><td>x/(exp(θ<sub>0</sub>) - |θ<sub>1</sub>|√x)</td><td>9</td><td>2.2 × 10<sup>-11</sup></td><td>0.03</td><td>0.44</td><td>-</td><td>-</td><td>-1250.9</td><td>17.5</td><td>7.3</td><td>-1226.1</td></tr><tr><td colspan="12" style="text-align: center;">1 - log L(θ̂) 2klog(n) + Σ<sub>j</sub>log(c<sub>j</sub>) 3 - p/2 log(3) + Σ<sub>i</sub>p log(I<sub>ii</sub>)<sup>1/2</sup> + log(|θ̂<sub>i</sub>|)</td></tr></tbody></table>

MOND cosmology and maximising agreement with the rotation curve fits; Chae et al. (2021), and hence providing an upper bound on the impact of the EFE. This is similar to the value inferred in Chae (2022) and Chae et al. (2022) from fits to the SPARC rotation curves, and from our fit to the SPARC data in Table 1.<fnref n="5" />

# 4 RESULTS

## 4.1 SPARC data

We show in Table 1 the statistics of the best functions found by ESR on the SPARC data. We split the codelength of equation (1) into terms describing the residuals of the data around the functional expectation, the functional form, and the parameter values as shown in the table footnotes. Below the horizontal line we give the results of the three MOND IFs, for which the free parameter corresponds to $g_0$, the double power law (equation (5)), and the Simple IF + EFE (equation (6)). $P(f) \equiv \exp(-L(D))/\sum(\exp(-L(D)))$ is the probability of the function given its description length, where the sum is over all functions up to complexity 9. (Note that these values would be changed by low-$L(D)$ functions at higher complexity.) For reference, $L(D)$ for the raw data (corresponding to the hypothesis $\log g_{obs} = 0$) is 53471, showing that significant compression is possible.

We find the best-fitting $g_0$ value for the RAR IF to be 1.13, somewhat lower than the 1.20 quoted by Lelli et al. (2017), although the data are the same. This is because Lelli et al. (2017) used

scipy.odr to perform the optimization rather than using the full first-order likelihood equations (2)–(3), and also did the analysis in the $\log(g_{obs}) - g_{bar}$ rather than $\log(g_{bar}) - \log(g_{obs})$ plane (F. Lelli private communication). The double power law fit has much higher maximum likelihood than the MOND IFs, outweighing its increased codelength due to its four free parameters. Although the RAR IF has complexity 9, it is not explicitly produced by ESR due to the constant ‘1’ appearing,<fnref n="6" /> a generalized form in which this is replaced by a free parameter appears at rank 17, with a probability $4 \times 10^{10}$ times lower than the top-ranked function ($\Delta L(D) = 24.5$). When $\theta_0 \neq 0$, the low-$g_{bar}$ logarithmic slope $s_-$ of this function is 1 rather than 1/2, so it does not function as a MOND IF. We refer to it as the ‘generalized RAR IF’.

The best ESR functions are clearly superior to the MOND functions or double power law. While the best metric for this is $L(D)$ (or equivalently $P(f)$), other statistics lead to the same conclusion. There are many functions more accurate (lower $-\log(L)$) than even the double power law. Although the functions at rank 1–9 have more free parameters than the IFs this is more than compensated for by their greater accuracy: as an alternative metric, the Bayesian Information Criterion (BIC) of the rank 1 function is 108 lower than the Simple IF, and even that of the rank 3 function with four free parameters is 94.5 lower, corresponding to a very strong preference. In agreement with Chae et al. (2020b), Chae et al. (2021), Chae et al. (2022) when fitting the Simple IF + EFE, we find that $e_N > 0$ is clearly preferred and recovers a value around the large-scale structure expectation. Although this function is significantly more accurate than any IF on its own, more than compensating for its additional free parameter

<footnote n="5">⁵Chae & Milgrom (2022) argue that equation (6) is only reliable with the inner points of galaxies’ rotation curves removed. As we are interested only in the approximate effect of the EFE on the low-$g_{bar}$ slope of the RAR, mainly sourced by outer rotation curve points, we do not apply a cut.</footnote>

<footnote n="6">⁶This will be changed in a future version of ESR that performs ‘integer snap’ of parameters where this reduces $L(D)$.</footnote>
<table><caption>**Table 1.** Top functions found by ESR applied to the SPARC data, ranked by total description length, compared to four MOND IFs, and a double power law below the horizontal line. $x \equiv g_{bar}/10^{-10} m s^{-2}$. We include also the 'generalized RAR IF' at rank 17, although this does not have a deep-MOND limit. The IFs and double power law are not produced explicitly by our implementation of ESR and hence their ranks are unknown, although they are clearly worse than the best ESR functions in both description length and likelihood. For the Simple, RAR and Standard IFs the parameter is $g_0$, while for the 'Simple IF + EFE' (equation 6), the first is $g_0$ and the second $e_N$.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th rowspan="2">Comp.</th><th rowspan="2">P(f)</th><th colspan="4">Parameters</th><th rowspan="2">Resid.<sup>1</sup></th><th colspan="2">Description length</th><th rowspan="2">Total</th></tr><tr><th>θ<sub>0</sub></th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th><th>θ<sub>3</sub></th><th>Func.<sup>2</sup></th><th>Param.<sup>3</sup></th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> (|θ<sub>1</sub> + x|<sup>2</sup> + x)</td><td>9</td><td>9.3 × 10<sup>-1</sup></td><td>0.84</td><td>-0.02</td><td>0.38</td><td>-</td><td>-1279.1</td><td>14.5</td><td>14.0</td><td>-1250.6</td></tr><tr><td>2</td><td>|θ<sub>1</sub>|<sup>x</sup> + θ<sub>0</sub>|<sup>2</sup> + x</td><td>9</td><td>6.4 × 10<sup>-2</sup></td><td>-0.99</td><td>0.64</td><td>0.36</td><td>-</td><td>-1279.9</td><td>12.5</td><td>19.6</td><td>-1247.9</td></tr><tr><td>3</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub> - x|<sup>2</sup> - θ<sub>3</sub></sup></td><td>9</td><td>2.0 × 10<sup>-3</sup></td><td>-1.4 × 10<sup>2</sup></td><td>0.02</td><td>0.14</td><td>0.89</td><td>-1276.4</td><td>12.5</td><td>19.5</td><td>-1244.4</td></tr><tr><td>4</td><td>|θ<sub>0</sub>(θ<sub>1</sub> + x)|<sup>2</sup> + x</td><td>9</td><td>1.4 × 10<sup>-4</sup></td><td>0.35</td><td>-0.02</td><td>0.34</td><td>-</td><td>-1268.9</td><td>14.5</td><td>12.7</td><td>-1241.7</td></tr><tr><td>5</td><td>|θ<sub>0</sub> - |θ<sub>1</sub> - x|<sup>2</sup>|<sup>3</sup></td><td>9</td><td>1.0 × 10<sup>-5</sup></td><td>-0.30</td><td>0.02</td><td>0.42</td><td>2.14</td><td>-1271.1</td><td>12.5</td><td>19.5</td><td>-1239.1</td></tr><tr><td>6</td><td>√x exp (|θ<sub>0</sub>+x|<sup>1</sup>/2)</td><td>9</td><td>1.5 × 10<sup>-9</sup></td><td>-0.02</td><td>0.36</td><td>-</td><td>-</td><td>-1257.9</td><td>17.5</td><td>10.0</td><td>-1230.3</td></tr><tr><td>7</td><td>(|θ<sub>0</sub>|<sup>x</sup>/x) + x</td><td>9</td><td>2.4 × 10<sup>-10</sup></td><td>1.87</td><td>-0.52</td><td>-</td><td>-</td><td>-1250.6</td><td>14.5</td><td>7.6</td><td>-1228.5</td></tr><tr><td>8</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>8</td><td>1.8 × 10<sup>-10</sup></td><td>-1.8 × 10<sup>-3</sup></td><td>0.72</td><td>-</td><td>-</td><td>-1245.6</td><td>12.9</td><td>4.5</td><td>-1228.2</td></tr><tr><td>9</td><td>|θ<sub>0</sub> + 1/√x|<sup>θ<sub>1</sub></sup></td><td>8</td><td>9.6 × 10<sup>-11</sup></td><td>-0.22</td><td>-2.14</td><td>-</td><td>-</td><td>-1251.1</td><td>14.3</td><td>9.2</td><td>-1227.6</td></tr><tr><td>10</td><td>(√x + 1/√x)<sup>θ<sub>0</sub></sup> + x</td><td>9</td><td>8.2 × 10<sup>-11</sup></td><td>-0.53</td><td>-</td><td>-</td><td>-</td><td>-1248.3</td><td>16.1</td><td>4.8</td><td>-1227.4</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>17</td><td>x/(exp(θ<sub>0</sub>) - |θ<sub>1</sub>|√x)</td><td>9</td><td>2.2 × 10<sup>-11</sup></td><td>0.03</td><td>0.44</td><td>-</td><td>-</td><td>-1250.9</td><td>17.5</td><td>7.3</td><td>-1226.1</td></tr><tr><td colspan="13"><sup>1</sup> - log L(θ̂) <sup>2</sup>klog(n) + Σ<sub>j</sub>log(c<sub>j</sub>) <sup>3</sup> - p/2 log(3) + Σ<sub>i</sub>p log(I<sub>ii</sub>)<sup>1/2</sup> + log(|θ̂<sub>i</sub>|)</td></tr></tbody></table>

MOND cosmology and maximising agreement with the rotation curve fits; Chae et al. (2021), and hence providing an upper bound on the impact of the EFE. This is similar to the value inferred in Chae (2022) and Chae et al. (2022) from fits to the SPARC rotation curves, and from our fit to the SPARC data in Table 1.<fnref n="5" />

# 4 RESULTS

## 4.1 SPARC data

We show in Table 1 the statistics of the best functions found by ESR on the SPARC data. We split the codelength of equation (1) into terms describing the residuals of the data around the functional expectation, the functional form, and the parameter values as shown in the table footnotes. Below the horizontal line we give the results of the three MOND IFs, for which the free parameter corresponds to $g_0$, the double power law (equation (5)), and the Simple IF + EFE (equation (6)). $P(f) \equiv \exp(-L(D))/\sum(\exp(-L(D)))$ is the probability of the function given its description length, where the sum is over all functions up to complexity 9. (Note that these values would be changed by low-$L(D)$ functions at higher complexity.) For reference, $L(D)$ for the raw data (corresponding to the hypothesis $\log g_{obs} = 0$) is 53471, showing that significant compression is possible.

We find the best-fitting $g_0$ value for the RAR IF to be 1.13, somewhat lower than the 1.20 quoted by Lelli et al. (2017), although the data are the same. This is because Lelli et al. (2017) used

scipy.odr to perform the optimization rather than using the full first-order likelihood equations (2)–(3), and also did the analysis in the $\log(g_{obs}) - g_{bar}$ rather than $\log(g_{bar}) - \log(g_{obs})$ plane (F. Lelli private communication). The double power law fit has much higher maximum likelihood than the MOND IFs, outweighing its increased codelength due to its four free parameters. Although the RAR IF has complexity 9, it is not explicitly produced by ESR due to the constant ‘1’ appearing,<fnref n="6" /> a generalized form in which this is replaced by a free parameter appears at rank 17, with a probability $4 \times 10^{10}$ times lower than the top-ranked function ($\Delta L(D) = 24.5$). When $\theta_0 \neq 0$, the low-$g_{bar}$ logarithmic slope $s_-$ of this function is 1 rather than 1/2, so it does not function as a MOND IF. We refer to it as the ‘generalized RAR IF’.

The best ESR functions are clearly superior to the MOND functions or double power law. While the best metric for this is $L(D)$ (or equivalently $P(f)$), other statistics lead to the same conclusion. There are many functions more accurate (lower $-\log(L)$) than even the double power law. Although the functions at rank 1–9 have more free parameters than the IFs this is more than compensated for by their greater accuracy: as an alternative metric, the Bayesian Information Criterion (BIC) of the rank 1 function is 108 lower than the Simple IF, and even that of the rank 3 function with four free parameters is 94.5 lower, corresponding to a very strong preference. In agreement with Chae et al. (2020b), Chae et al. (2021), Chae et al. (2022) when fitting the Simple IF + EFE, we find that $e_N > 0$ is clearly preferred and recovers a value around the large-scale structure expectation. Although this function is significantly more accurate than any IF on its own, more than compensating for its additional free parameter

<footnote n="5">⁵Chae & Milgrom (2022) argue that equation (6) is only reliable with the inner points of galaxies’ rotation curves removed. As we are interested only in the approximate effect of the EFE on the low-$g_{bar}$ slope of the RAR, mainly sourced by outer rotation curve points, we do not apply a cut.</footnote>

<footnote n="6">⁶This will be changed in a future version of ESR that performs ‘integer snap’ of parameters where this reduces $L(D)$.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 2.** The top 20 functions found by ESR overlaid on the SPARC data (blue points), colour-coded by their relative probability in the full function list. The top panel fits the SPARC data, the middle panel mock data generated by the RAR IF, and the bottom panel mock data generated by the Simple IF with universal external field strength $g_{ex} = 1.2 \times 10^{-12} \text{ m s}^{-2}$. The mock data sets are 10 times larger than SPARC, although this is factored out in the description length calculation.</figcaption></figure>

($\Delta$BIC = -35 compared to the Simple IF), it has a poor description length due to the large functional contribution. The complexity is 59 using our current basis set of operators, although this would fall to 45 if tanh were explicitly included. In general, the great improvement in accuracy and simplicity of the ESR functions demonstrates the advantage of this method over guessing functions 'by eye'.

In the top panel of Fig. 2, we plot the best 20 functions on top of the SPARC data. The functions are colour-coded by $P(f)$, with darker colouring indicating functions favoured by MDL. The top six functions have discontinuities in $s$ at $g_{bar} \approx 0.02$. We have checked that this is not due to outlying points, does not invalidate the first-order approximation of equation (2), and does not appear to map onto any local feature of the data. Instead it is likely due to the relative simplicity (i.e. complexity $\le 9$) of the functions considered, as we discuss further in Section 5. While such behaviour could be excluded by requiring no discontinuity in the derivative within the range of the data (or more generally weighting functions with an $s$-dependent prior), we see no principled reason to do so. The first function without a discontinuity is at rank 7, which has the Newtonian limit and $s_- = 0.52$ at maximum likelihood. Although this is $9.3 \sigma$ from $s_- = 1/2$, and hence this equation cannot function as a pure-MOND IF, the

EFE leads to the expectation that $s > 1/2$ at low $g_{bar}$ as discussed in more detail below. Several of the remaining highly ranked functions have similar quasi-MONDian behaviour.

While some of the best functions have MONDian limits around $g_{bar} = 10^{\pm 3}$ others do not, and in particular $s_- < 1/2$ is common. To explore this further, we calculate $s_-$ and $s_+$ analytically for the top 10 equations as a function of their free parameters, showing the results in Table B1. For each of the functions where $s_-$ and/or $s_+$ depend on $\theta$ (i.e. are not fixed by the functional form alone), we perform a Markov Chain Monte Carlo (MCMC) inference using the `numpyro` sampler (Bingham et al. 2019; Phan, Pradhan & Jankowiak 2019) with broad flat priors to constrain the parameters, and hence derive the posterior predictive distributions of the limiting slopes. Fig. 3 (left-hand panel) shows the results for the top 10 functions, using a dot to indicate a limit fixed by the functional form, a bar to show the 95 per cent confidence interval in cases where the slope depends on the parameters, and an arrowhead to indicate a limit outside of the plotting range (including $\pm\infty$). In many cases, the 95 per cent confidence interval is extremely narrow.

Although it is $s_-$ and $s_+$ that directly relate to the MOND hypothesis, they require extrapolation far beyond the range of the data. To understand how the slopes behave near the limits of the SPARC data, we also calculate $s$ at the minimum, $g_{bar, min} = 8.32 \times 10^{-13} \text{ m s}^{-2}$, and maximum, $g_{bar, max} = 6.54 \times 10^{-9} \text{ m s}^{-2}$, measured baryonic accelerations. These are plotted in cyan and magenta respectively in Fig. 3. At $g_{bar, max}$, the logarithmic slope is $\lesssim 1$ for almost all functions, as expected for the Newtonian limit as only at $g_{bar} \to \infty$ does $s$ become 1. However, we find that the $g_{bar, min}$ slopes of the top five functions are not $\sim 1/2$ but actually $< 0$ due to the aforementioned discontinuity. The remaining top functions have low-acceleration slopes typically slightly larger than $1/2$.

These results show that the SPARC data do not unambiguously favour $s_- = 1/2$ and $s_+ = 1$. The requirement for an interpolating function to be MONDian is in fact even more stringent than this, since the coefficient of proportionality in the limiting high-$g_{bar}$ power-law relation must be unity, i.e. $g_{obs} = g_{bar}$. We find that among the functions in the top 10 for which $s_+ = 1$, four have such a coefficient (at rank 2, 4, 7, 10) while for the rank 1 function this is $0.84 \pm 0.006$ and at rank 8 it is $0.72 \pm 0.01$, where the uncertainties are obtained by fitting the functions with MCMC. At low $g_{bar}$, the coefficient in $g_{obs} \propto g_{bar}^{1/2}$ should be $\sqrt{g_0} \approx 1$. The only function with $s_- = 1/2$ (rank 6) has the coefficient $1.12$ (close to the 1.10 expected from the canonical $g_0 = 1.2$) while those further down with $s_{gbar, min} \approx 1/2$ have a coefficient of 1. The relative simplicity of the functions we consider here should preference a coefficient of 1, and hence again it is not clear to what extent the data may be said to be MONDian. The double power law has the limits $g_{obs} = 0.81 g_{bar}^{1.03}$ at high $g_{bar}$ and $g_{obs} = 1.57 g_{bar}^{0.60}$ at low $g_{bar}$.

Next, we show in the left-hand panel of Fig. 4, the separate Pareto fronts of description length and negative log-likelihood, with the second ('simplicity') objective measured by functional complexity. Unlike the Pareto fronts produced by traditional SR algorithms, those of ESR are guaranteed to be optimal. $L(D)$ and $-\log(L)$ are minimized separately at each complexity, and have their minimum values over all complexity subtracted so that the globally best functions appear at 0. We show the MONDian functions and double power law as separate symbols, all of which we find to be strongly Pareto-dominated by the best ESR functions at lower complexity. Note that while the 'knee' of the Pareto front (where $L(D)$ or $-\log(L)$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** The top 20 functions found by ESR overlaid on the SPARC data (blue points), colour-coded by their relative probability in the full function list. The top panel fits the SPARC data, the middle panel mock data generated by the RAR IF, and the bottom panel mock data generated by the Simple IF with universal external field strength $g_{ex} = 1.2 \times 10^{-12} \text{ m s}^{-2}$. The mock data sets are 10 times larger than SPARC, although this is factored out in the description length calculation.</figcaption></figure>

($\Delta$BIC = -35 compared to the Simple IF), it has a poor description length due to the large functional contribution. The complexity is 59 using our current basis set of operators, although this would fall to 45 if tanh were explicitly included. In general, the great improvement in accuracy and simplicity of the ESR functions demonstrates the advantage of this method over guessing functions 'by eye'.

In the top panel of Fig. 2, we plot the best 20 functions on top of the SPARC data. The functions are colour-coded by $P(f)$, with darker colouring indicating functions favoured by MDL. The top six functions have discontinuities in $s$ at $g_{bar} \approx 0.02$. We have checked that this is not due to outlying points, does not invalidate the first-order approximation of equation (2), and does not appear to map onto any local feature of the data. Instead it is likely due to the relative simplicity (i.e. complexity $\le 9$) of the functions considered, as we discuss further in Section 5. While such behaviour could be excluded by requiring no discontinuity in the derivative within the range of the data (or more generally weighting functions with an $s$-dependent prior), we see no principled reason to do so. The first function without a discontinuity is at rank 7, which has the Newtonian limit and $s_- = 0.52$ at maximum likelihood. Although this is $9.3 \sigma$ from $s_- = 1/2$, and hence this equation cannot function as a pure-MOND IF, the

EFE leads to the expectation that $s > 1/2$ at low $g_{bar}$ as discussed in more detail below. Several of the remaining highly ranked functions have similar quasi-MONDian behaviour.

While some of the best functions have MONDian limits around $g_{bar} = 10^{\pm 3}$ others do not, and in particular $s_- < 1/2$ is common. To explore this further, we calculate $s_-$ and $s_+$ analytically for the top 10 equations as a function of their free parameters, showing the results in Table B1. For each of the functions where $s_-$ and/or $s_+$ depend on $\theta$ (i.e. are not fixed by the functional form alone), we perform a Markov Chain Monte Carlo (MCMC) inference using the `numpyro` sampler (Bingham et al. 2019; Phan, Pradhan & Jankowiak 2019) with broad flat priors to constrain the parameters, and hence derive the posterior predictive distributions of the limiting slopes. Fig. 3 (left-hand panel) shows the results for the top 10 functions, using a dot to indicate a limit fixed by the functional form, a bar to show the 95 per cent confidence interval in cases where the slope depends on the parameters, and an arrowhead to indicate a limit outside of the plotting range (including $\pm\infty$). In many cases, the 95 per cent confidence interval is extremely narrow.

Although it is $s_-$ and $s_+$ that directly relate to the MOND hypothesis, they require extrapolation far beyond the range of the data. To understand how the slopes behave near the limits of the SPARC data, we also calculate $s$ at the minimum, $g_{bar, min} = 8.32 \times 10^{-13} \text{ m s}^{-2}$, and maximum, $g_{bar, max} = 6.54 \times 10^{-9} \text{ m s}^{-2}$, measured baryonic accelerations. These are plotted in cyan and magenta respectively in Fig. 3. At $g_{bar, max}$, the logarithmic slope is $\lesssim 1$ for almost all functions, as expected for the Newtonian limit as only at $g_{bar} \to \infty$ does $s$ become 1. However, we find that the $g_{bar, min}$ slopes of the top five functions are not $\sim 1/2$ but actually $< 0$ due to the aforementioned discontinuity. The remaining top functions have low-acceleration slopes typically slightly larger than $1/2$.

These results show that the SPARC data do not unambiguously favour $s_- = 1/2$ and $s_+ = 1$. The requirement for an interpolating function to be MONDian is in fact even more stringent than this, since the coefficient of proportionality in the limiting high-$g_{bar}$ power-law relation must be unity, i.e. $g_{obs} = g_{bar}$. We find that among the functions in the top 10 for which $s_+ = 1$, four have such a coefficient (at rank 2, 4, 7, 10) while for the rank 1 function this is $0.84 \pm 0.006$ and at rank 8 it is $0.72 \pm 0.01$, where the uncertainties are obtained by fitting the functions with MCMC. At low $g_{bar}$, the coefficient in $g_{obs} \propto g_{bar}^{1/2}$ should be $\sqrt{g_0} \approx 1$. The only function with $s_- = 1/2$ (rank 6) has the coefficient $1.12$ (close to the 1.10 expected from the canonical $g_0 = 1.2$) while those further down with $s_{gbar, min} \approx 1/2$ have a coefficient of 1. The relative simplicity of the functions we consider here should preference a coefficient of 1, and hence again it is not clear to what extent the data may be said to be MONDian. The double power law has the limits $g_{obs} = 0.81 g_{bar}^{1.03}$ at high $g_{bar}$ and $g_{obs} = 1.57 g_{bar}^{0.60}$ at low $g_{bar}$.

Next, we show in the left-hand panel of Fig. 4, the separate Pareto fronts of description length and negative log-likelihood, with the second ('simplicity') objective measured by functional complexity. Unlike the Pareto fronts produced by traditional SR algorithms, those of ESR are guaranteed to be optimal. $L(D)$ and $-\log(L)$ are minimized separately at each complexity, and have their minimum values over all complexity subtracted so that the globally best functions appear at 0. We show the MONDian functions and double power law as separate symbols, all of which we find to be strongly Pareto-dominated by the best ESR functions at lower complexity. Note that while the 'knee' of the Pareto front (where $L(D)$ or $-\log(L)$
<figure><img src="image_3.png" /><figcaption>**Figure 3.** The logarithmic slopes $s \equiv \frac{d\log(g_{\text{obs}})}{d\log(g_{\text{bar}})}$ of the top 10 ESR functions on each data set, for comparison with the low- and high-$g_{\text{bar}}$ MONDian expectations 1/2 and 1, respectively (blue and red vertical dashed lines). The blue and red points are the limiting slopes $s_{-} \equiv \lim_{g_{\text{bar}} \to 0^{+}} s$ and $s_{+} \equiv \lim_{g_{\text{bar}} \to \infty} s$, while cyan and magenta indicate the slopes at the minimum and maximum $g_{\text{bar}}$ of the SPARC data (0.0083 and 65.4). In case a slope depends on a parameter value, we show the 95 per cent confidence interval as a bar (often very thin) obtained from an MCMC fit. Arrowheads indicate points or bars beyond the range of the plot.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 4.** The Pareto fronts identified by ESR for the SPARC, RAR IF mock, and Simple IF + EFE mock data sets, for both $\log(\mathcal{L})$ (blue) and total description length $L(D)$ (red). The quantities plotted have the minimum values subtracted so that the best results appear at 0. Also shown are the results of the RAR, Simple and Standard IFs, Simple IF + EFE, and double power law fits. ESR significantly outperforms these ‘by eye’ guesses, even for mock data generated from them. Short diagonal lines on the x-axis indicate breaks. In the left-hand and right-hand panels both red and blue points for the Standard IF at complexity 14 lie above the top of the plot.</figcaption></figure>

turns over) would appear to be at complexity 6–7, there is a significant improvement in going from complexity 8 to 9. This cautions against automatically selecting functions at the knee (the default for example in PYSR), and indicates that further improvement would likely be achievable by going beyond complexity 9. This is beyond the scope of the present work; we are content here to have discovered simple functional forms for the RAR surpassing any that have been considered heretofore.

## 4.2 Mock data

The above results suggest a Newtonian limit ($g_{\text{obs}} \to g_{\text{bar}}$ as $g_{\text{bar}} \to \infty$) is somewhat favoured by the data while a deep-MOND limit ($g_{\text{obs}} \to \sqrt{g_0} g_{\text{bar}}$ as $g_{\text{bar}} \to 0$) is questionable. However, given the

limited dynamical range and significant uncertainties of the data, it is unclear to what extent we should expect to find these limits even if the generating function were MONDian. In addition, the EFE would imply $s > 1/2$ at low $g_{\text{bar}}$. To investigate these issues we now apply ESR to the mock data of Section 3.4.

### 4.2.1 RAR IF generating function

Table 2 shows the best functions found by ESR for the RAR IF mock data, along with the results for the IFs and double power law. The RAR IF is, by construction, a good fit to this data, but there are 15 functions with lower $L(D)$, including several at lower complexity. This indicates that the characteristics of the SPARC data (dynamic range and uncertainties) are insufficient to pick out
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The logarithmic slopes $s \equiv \frac{d\log(g_{\text{obs}})}{d\log(g_{\text{bar}})}$ of the top 10 ESR functions on each data set, for comparison with the low- and high-$g_{\text{bar}}$ MONDian expectations 1/2 and 1, respectively (blue and red vertical dashed lines). The blue and red points are the limiting slopes $s_{-} \equiv \lim_{g_{\text{bar}} \to 0^{+}} s$ and $s_{+} \equiv \lim_{g_{\text{bar}} \to \infty} s$, while cyan and magenta indicate the slopes at the minimum and maximum $g_{\text{bar}}$ of the SPARC data (0.0083 and 65.4). In case a slope depends on a parameter value, we show the 95 per cent confidence interval as a bar (often very thin) obtained from an MCMC fit. Arrowheads indicate points or bars beyond the range of the plot.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** The Pareto fronts identified by ESR for the SPARC, RAR IF mock, and Simple IF + EFE mock data sets, for both $\log(\mathcal{L})$ (blue) and total description length $L(D)$ (red). The quantities plotted have the minimum values subtracted so that the best results appear at 0. Also shown are the results of the RAR, Simple and Standard IFs, Simple IF + EFE, and double power law fits. ESR significantly outperforms these ‘by eye’ guesses, even for mock data generated from them. Short diagonal lines on the x-axis indicate breaks. In the left-hand and right-hand panels both red and blue points for the Standard IF at complexity 14 lie above the top of the plot.</figcaption></figure>

turns over) would appear to be at complexity 6–7, there is a significant improvement in going from complexity 8 to 9. This cautions against automatically selecting functions at the knee (the default for example in PYSR), and indicates that further improvement would likely be achievable by going beyond complexity 9. This is beyond the scope of the present work; we are content here to have discovered simple functional forms for the RAR surpassing any that have been considered heretofore.

## 4.2 Mock data

The above results suggest a Newtonian limit ($g_{\text{obs}} \to g_{\text{bar}}$ as $g_{\text{bar}} \to \infty$) is somewhat favoured by the data while a deep-MOND limit ($g_{\text{obs}} \to \sqrt{g_0} g_{\text{bar}}$ as $g_{\text{bar}} \to 0$) is questionable. However, given the

limited dynamical range and significant uncertainties of the data, it is unclear to what extent we should expect to find these limits even if the generating function were MONDian. In addition, the EFE would imply $s > 1/2$ at low $g_{\text{bar}}$. To investigate these issues we now apply ESR to the mock data of Section 3.4.

### 4.2.1 RAR IF generating function

Table 2 shows the best functions found by ESR for the RAR IF mock data, along with the results for the IFs and double power law. The RAR IF is, by construction, a good fit to this data, but there are 15 functions with lower $L(D)$, including several at lower complexity. This indicates that the characteristics of the SPARC data (dynamic range and uncertainties) are insufficient to pick out
The middle panel of Fig. 4 shows the Pareto front for these data. We find more smooth behaviour than for the real data, with the optimum solution achieved already by complexity 8, reflecting the relatively simple nature of the generating function. This shows that if the RAR IF was generating the real data (and our likelihood and mock data generation method were accurate), we would have achieved the $L(D)$ minimum on those data too. However, the MONDian functions (including the RAR IF itself) and double power law are Pareto-dominated by the ESR results even on these mock data, showing that one should not expect to be able to recover unambiguously even this simplest of MOND generating functions.

### 4.2.2 Simple IF + EFE generating function

Analogous results for the mock data generated using the Simple IF with inclusion of the EFE are shown in Table 3 and the bottom/right-hand panels of Figs 2–4. This data set behaves more similarly to the real data in terms of the relative ordering of the IFs, double power law, and ESR functions, including the generalized RAR IF. Indeed, the best-fit parameters of all three non-EFE IFs are identical to the SPARC data to two decimal places, while those of the generalized RAR IF and the high-$g_{bar}$ slope of the double power law are the same to one. Here, the IFs provide a significantly worse compression of the data than the best ESR functions, and the double power law also performs relatively poorly due to the curvature at low $g_{bar}$ (see Fig. 1). There is again a small bias between the maximum-likelihood (1.19 and $8.56 \times 10^{-3}$) and true (1.2 and $1 \times 10^{-2}$) $g_0$ and $e_N$ values for the Simple IF + EFE fit. Although this function has among the highest likelihoods achievable by ESR up to complexity 9, its functional complexity makes it a poor compression of its own SPARC-like mock data. This reinforces the conclusion that the characteristics of these data are insufficient to identify a MONDian generating function: this one in particular would require far more data than the RAR IF to be favoured by MDL.

For the Simple IF + EFE mock all top-10 functions have $s_+ = 1$ and $s(g_{bar,max}) > 0.9$. However, only six of them recover the true Newtonian limit $g_{obs} = g_{bar}$; the others find $g_{obs} = 0.71 g_{bar}$ or $g_{obs} = 0.79 g_{bar}$, again with sub percent uncertainty from MCMC. Thus under this model too one would not expect the Newtonian limit to be identified robustly. While $s(g_{bar,min}) > 1/2$, indicating the significant impact of the EFE, $s_-$ is typically 0 as opposed to 1 as expected from equation (6). Thus $g_{bar,min}$ is too high to constrain $s_-$ reliably, although this may also be a reflection of the relative simplicity of the functions we consider. The double power law limits are $g_{obs} = 0.96 g_{bar}^{0.98}$ at high $g_{bar}$ and $g_{obs} = 1.55 g_{bar}^{0.60}$ at low $g_{bar}$. The Pareto front indicates this data set to be somewhat more complex than the RAR IF mock, as $L(D)$ continues to fall to complexity 9, although the smoothness shows it to be simpler than the real data.

All functions considered achieve considerably higher likelihood on the mock data sets than the real data, showing that the mocks are simpler. This could be because the data do not conform to the MOND expectation – one would expect any given $g_{obs} = f(g_{bar})$ to be relatively inaccurate in a chaotic ΛCDM galaxy formation scenario – or because the model for scattering the mock data points is overly simplistic. This is discussed further in Section 5. Relatedly, the $P(f)$ values of the top functions are very closely spaced in the mock data sets, indicating that there is little to distinguish them. On the contrary, on the real data the integrated probability of all functions besides the top five is $\lesssim 10^{-8}$, suggesting that these functions perhaps ought not to be considered at all. The limiting slopes of the top 10 equations as

a function of the parameters can be found in Tables B2 and B3 for the RAR and EFE mocks, respectively.

## 5 DISCUSSION

Our main conclusion is that the SPARC data are insufficient to determine robustly the limiting behaviour of the RAR, and hence cannot verify or refute the MOND hypothesis. This is reached by studying mock data generated by MOND; in particular, generating data according to the RAR IF, not only are we unable to identify that as the generating function, but, more seriously, we cannot reconstruct $s_- = 1/2$. At the high-$g_{bar}$ end, the logarithmic slope of the Newtonian limit ($s_+ = 1$) is typically well recovered, although the coefficient of proportionality in $g_{obs} \propto g_{bar}$ is not: in the RAR IF mock data this takes a values $\sim 0.64$ far more often than 1.

Improving this situation requires increasing the dynamical range of the RAR. At the low-$g_{bar}$ end, this may be achieved by studying ultra-diffuse galaxies (e.g. Freundlich et al. 2022), or local dwarf spheroidals (e.g. McGaugh & Wolf 2010; McGaugh & Milgrom 2013), some of which seem to indicate $s_- \approx 0$, as found by many of the best ESR functions (Lelli et al. 2017). Alternatively, one may attempt to probe the outer regions of galaxies including the Milky Way (e.g. Oman et al. 2020). Particularly promising for a large gain is to use stacked weak lensing to probe galaxy outskirts that would have insufficient signal to noise on an individual-object basis (Brouwer et al. 2021). This appears to indicate $s_- \approx 1/2$. Increasing $g_{bar,max}$ requires probing the central regions of high-mass ellipticals, as well as groups and clusters of galaxies (Chae et al. 2019, 2020a; Chan & Del Popolo 2020; Tian et al. 2020; Gopika & Desai 2021; Pradyumna & Desai 2021). Such data already exists and may readily be folded into our framework to increase its constraining power. A smaller information gain may be achieved by reducing the uncertainties in $g_{bar}$ and $g_{obs}$: in the limit of no uncertainty, any generating function will be assigned $P(f) = 1$ by MDL. By generating mock data with different characteristics, one could ascertain the requirements for various features of the functional form to be unambiguously determined.

It is likely that there exist functions at higher complexity superior to those of Tables 1–3, especially for the real data where $L(D)$ drops significantly from complexity 8 to 9. While uncovering lower-$L(D)$ functions at complexity >9 may update the optimal limiting behaviours of the functional form of the RAR, and hence its compatibility with MOND, it cannot compromise our discovery of simpler and more accurate functions than the IFs and double power law. Indeed, the fact that the (or at least *a*) knee of the Pareto front is reached around complexity 7 in the left-hand panel of Fig. 4 shows that such functions already offer a powerful compression, and the commonalities between the top functions at complexity $\lesssim 9$ suggests that similar features are likely to be present in more complex functions also.

Discovering such functions would likely be computationally prohibitive for the ESR algorithm, and thus a stochastic search (e.g. using a genetic algorithm) may be required. This search may be seeded by the ESR functions: the fact that many of the best-fitting functions have similar features (such as $g_{obs} \approx 0.64 g_{bar}$ as $g_{bar} \to \infty$ for the RAR IF mock) suggests these may be useful for higher-complexity functions also. Thus ESR may be used to validate the underlying assumption of stochastic searches that there exist features of functions responsible for their fitness, and the identification of these features may be useful for tuning hyperparameters. It would also be possible to combine ESR with deterministic symbolic regression algorithms (e.g. Worm & Chiu 2013; Rivero & Fernandez-Blanco 2022; Kammerer et al.
The middle panel of Fig. 4 shows the Pareto front for these data. We find more smooth behaviour than for the real data, with the optimum solution achieved already by complexity 8, reflecting the relatively simple nature of the generating function. This shows that if the RAR IF was generating the real data (and our likelihood and mock data generation method were accurate), we would have achieved the $L(D)$ minimum on those data too. However, the MONDian functions (including the RAR IF itself) and double power law are Pareto-dominated by the ESR results even on these mock data, showing that one should not expect to be able to recover unambiguously even this simplest of MOND generating functions.

### 4.2.2 Simple IF + EFE generating function

Analogous results for the mock data generated using the Simple IF with inclusion of the EFE are shown in Table 3 and the bottom/right-hand panels of Figs 2–4. This data set behaves more similarly to the real data in terms of the relative ordering of the IFs, double power law, and ESR functions, including the generalized RAR IF. Indeed, the best-fit parameters of all three non-EFE IFs are identical to the SPARC data to two decimal places, while those of the generalized RAR IF and the high-$g_{bar}$ slope of the double power law are the same to one. Here, the IFs provide a significantly worse compression of the data than the best ESR functions, and the double power law also performs relatively poorly due to the curvature at low $g_{bar}$ (see Fig. 1). There is again a small bias between the maximum-likelihood (1.19 and $8.56 \times 10^{-3}$) and true (1.2 and $1 \times 10^{-2}$) $g_0$ and $e_N$ values for the Simple IF + EFE fit. Although this function has among the highest likelihoods achievable by ESR up to complexity 9, its functional complexity makes it a poor compression of its own SPARC-like mock data. This reinforces the conclusion that the characteristics of these data are insufficient to identify a MONDian generating function: this one in particular would require far more data than the RAR IF to be favoured by MDL.

For the Simple IF + EFE mock all top-10 functions have $s_+ = 1$ and $s(g_{bar,max}) > 0.9$. However, only six of them recover the true Newtonian limit $g_{obs} = g_{bar}$; the others find $g_{obs} = 0.71 g_{bar}$ or $g_{obs} = 0.79 g_{bar}$, again with sub percent uncertainty from MCMC. Thus under this model too one would not expect the Newtonian limit to be identified robustly. While $s(g_{bar,min}) > 1/2$, indicating the significant impact of the EFE, $s_-$ is typically 0 as opposed to 1 as expected from equation (6). Thus $g_{bar,min}$ is too high to constrain $s_-$ reliably, although this may also be a reflection of the relative simplicity of the functions we consider. The double power law limits are $g_{obs} = 0.96 g_{bar}^{0.98}$ at high $g_{bar}$ and $g_{obs} = 1.55 g_{bar}^{0.60}$ at low $g_{bar}$. The Pareto front indicates this data set to be somewhat more complex than the RAR IF mock, as $L(D)$ continues to fall to complexity 9, although the smoothness shows it to be simpler than the real data.

All functions considered achieve considerably higher likelihood on the mock data sets than the real data, showing that the mocks are simpler. This could be because the data do not conform to the MOND expectation – one would expect any given $g_{obs} = f(g_{bar})$ to be relatively inaccurate in a chaotic ΛCDM galaxy formation scenario – or because the model for scattering the mock data points is overly simplistic. This is discussed further in Section 5. Relatedly, the $P(f)$ values of the top functions are very closely spaced in the mock data sets, indicating that there is little to distinguish them. On the contrary, on the real data the integrated probability of all functions besides the top five is $\lesssim 10^{-8}$, suggesting that these functions perhaps ought not to be considered at all. The limiting slopes of the top 10 equations as

a function of the parameters can be found in Tables B2 and B3 for the RAR and EFE mocks, respectively.

## 5 DISCUSSION

Our main conclusion is that the SPARC data are insufficient to determine robustly the limiting behaviour of the RAR, and hence cannot verify or refute the MOND hypothesis. This is reached by studying mock data generated by MOND; in particular, generating data according to the RAR IF, not only are we unable to identify that as the generating function, but, more seriously, we cannot reconstruct $s_- = 1/2$. At the high-$g_{bar}$ end, the logarithmic slope of the Newtonian limit ($s_+ = 1$) is typically well recovered, although the coefficient of proportionality in $g_{obs} \propto g_{bar}$ is not: in the RAR IF mock data this takes a values $\sim 0.64$ far more often than 1.

Improving this situation requires increasing the dynamical range of the RAR. At the low-$g_{bar}$ end, this may be achieved by studying ultra-diffuse galaxies (e.g. Freundlich et al. 2022), or local dwarf spheroidals (e.g. McGaugh & Wolf 2010; McGaugh & Milgrom 2013), some of which seem to indicate $s_- \approx 0$, as found by many of the best ESR functions (Lelli et al. 2017). Alternatively, one may attempt to probe the outer regions of galaxies including the Milky Way (e.g. Oman et al. 2020). Particularly promising for a large gain is to use stacked weak lensing to probe galaxy outskirts that would have insufficient signal to noise on an individual-object basis (Brouwer et al. 2021). This appears to indicate $s_- \approx 1/2$. Increasing $g_{bar,max}$ requires probing the central regions of high-mass ellipticals, as well as groups and clusters of galaxies (Chae et al. 2019, 2020a; Chan & Del Popolo 2020; Tian et al. 2020; Gopika & Desai 2021; Pradyumna & Desai 2021). Such data already exists and may readily be folded into our framework to increase its constraining power. A smaller information gain may be achieved by reducing the uncertainties in $g_{bar}$ and $g_{obs}$: in the limit of no uncertainty, any generating function will be assigned $P(f) = 1$ by MDL. By generating mock data with different characteristics, one could ascertain the requirements for various features of the functional form to be unambiguously determined.

It is likely that there exist functions at higher complexity superior to those of Tables 1–3, especially for the real data where $L(D)$ drops significantly from complexity 8 to 9. While uncovering lower-$L(D)$ functions at complexity $>9$ may update the optimal limiting behaviours of the functional form of the RAR, and hence its compatibility with MOND, it cannot compromise our discovery of simpler and more accurate functions than the IFs and double power law. Indeed, the fact that the (or at least *a*) knee of the Pareto front is reached around complexity 7 in the left-hand panel of Fig. 4 shows that such functions already offer a powerful compression, and the commonalities between the top functions at complexity $\lesssim 9$ suggests that similar features are likely to be present in more complex functions also.

Discovering such functions would likely be computationally prohibitive for the ESR algorithm, and thus a stochastic search (e.g. using a genetic algorithm) may be required. This search may be seeded by the ESR functions: the fact that many of the best-fitting functions have similar features (such as $g_{obs} \approx 0.64 g_{bar}$ as $g_{bar} \to \infty$ for the RAR IF mock) suggests these may be useful for higher-complexity functions also. Thus ESR may be used to validate the underlying assumption of stochastic searches that there exist features of functions responsible for their fitness, and the identification of these features may be useful for tuning hyperparameters. It would also be possible to combine ESR with deterministic symbolic regression algorithms (e.g. Worm & Chiu 2013; Rivero & Fernandez-Blanco 2022; Kammerer et al.
not to say either that data could not plausibly have been generated by MOND or that it could plausibly have been generated under another hypothesis, as only MOND currently has sufficient predictivity for a test of this precision. We look to future SR studies with more data to establish the functional form of the RAR – if it exists – definitively.

# 6 CONCLUSION

The RAR has become central to debates about the mass discrepancy problem on astrophysical scales. Its tightness and regularity have been used to argue for a violation of Newtonian gravity in accordance with MOND, but the functions used to fit the data have been constructed to conform to this theory. As the first detailed application of the brand-new technique of ESR, we rank objectively *all* simple functions in terms of their aptitude for describing the SPARC RAR. We employ the minimum description length principle to trade accuracy with simplicity and hence perform model selection, and calibrate our method on mock MOND data generated both with and without the external field effect (EFE). Our conclusions are as follows:

(i) ESR discovers functions which are better descriptions, in both accuracy and simplicity and for both observed and simulated data, than MOND functions or a double power law.

(ii) While the majority of best-fitting functions on the SPARC data recover $g_{\text{obs}} \propto g_{\text{bar}}$ at high accelerations, not all have a best-fit coefficient of proportionality near unity. Thus the Newtonian limit is not clearly evidenced.

(iii) The SPARC data do not prefer functions with the deep-MOND limit of $g_{\text{obs}} \propto \sqrt{g_{\text{bar}}}$ as $g_{\text{bar}} \to 0$. Instead, we find that functions with $g_{\text{bar}} \to \text{const}$ typically compress the data more efficiently, albeit with considerable uncertainty.

(iv) SPARC-like mock data generated assuming the MONDian RAR interpolating function do not unambiguously recover that function. Moreover, many of the best functions for those mock data have $g_{\text{obs}} \approx 0.64 g_{\text{bar}}$ rather than $g_{\text{obs}} = g_{\text{bar}}$ at high $g_{\text{bar}}$, and most do not have a deep-MOND limit at all.

(v) The EFE in AQUAL greatly increases the logarithmic slope of the best-fitting functions at the low-$g_{\text{bar}}$ end of the data, but does not appreciably impact the limiting slope at $g_{\text{bar}} \to 0$. Incorporating the EFE in the mock data produces more generally similar results to the real data, so our analysis (within the MOND paradigm) hints at it.

(vi) We conclude that the data have too small a dynamic range (and too large uncertainties) to unambiguously favour MOND even if it is in fact generating the data. The SPARC RAR alone therefore does not supporting that theory unambiguously. The best prospect for improving this situation is to increase the acceleration range of the data, e.g. using stacked weak lensing at low $g_{\text{bar}}$ and groups and clusters at high $g_{\text{bar}}$.

(vii) Our results are a function of the maximum complexity of equation considered. Future symbolic regression algorithms – exhaustive or non-exhaustive – will reach the true description length minimum, and hence uncover the optimal functional representation of the RAR and determine whether the relation implies novel law-like gravitational behaviour.

ESR provides for the first time a guaranteed complete search through functional parameter space, making it the ideal tool to determine the analytic form of observed relations, extract physics from data theory-agnostically, and create fitting functions. We make the ESR and RAR codes, full function sets and the best 50 functions for each data set we consider publicly available to facilitate future applications.

# ACKNOWLEDGEMENTS

We thank Kyu-Hyun Chae, Andrei Constantin, Miles Cranmer, Mario Figueiredo, Gianluca Gregori, Thomas Harvey, Mark Kotanchek, Federico Lelli, Stacy McGaugh, Andre Lukas, Richard Stiskalek, and Tariq Yasin for useful inputs and discussion.

HD is supported by a Royal Society University Research Fellowship (grant no. 211046). DJB is supported by the Simons Collaboration on ‘Learning the Universe’ and was supported by STFC and Oriel College Oxford. PGF acknowledges support from European Research Council Grant No: 693024 and the Beecroft Trust.

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

This work used the DiRAC Complexity and DiRAC@Durham facilities, operated by the University of Leicester IT Services and Institute for Computational Cosmology, which form part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by Business Innovation & Skills (BIS) National E-Infrastructure capital grants ST/K000373/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, STFC DiRAC Operations grant ST/K0003259/1, and Durham University and Science and Technology Facilities Council (STFC) operations grant ST/R000832/1. DiRAC is part of the National E-Infrastructure.

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

# 7 DATA AVAILABILITY

The code and data associated with ESR and its application to the RAR are released at https://github.com/DeaglanBartlett/ESR and in Bartlett, Desmond & Ferreira (2022a). The SPARC data is available at http://astroweb.cwru.edu/SPARC. Other data may be shared on request to the corresponding authors.

# REFERENCES

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

Bartlett D. J., Desmond H., Ferreira P. G., 2022a, Exhaustive Symbolic Regression Function Sets

Bartlett D. J., Desmond H., Ferreira P. G., 2022b, preprint (arXiv:2211.11461)

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

Berezhiani L., Khoury J., 2015, Phys. Rev. D, 92, 103510

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

Blanchet L., Le Tiec A., 2008, Phys. Rev. D, 78, 024031

Brouwer M. M. et al., 2021, A&A, 650, A113

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

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

Chae K.-H., Bernardi M., Sheth R. K., Gong I.-T., 2019, ApJ, 877, 18

Chae K.-H., Bernardi M., Domínguez Sánchez H., Sheth R. K., 2020a, ApJ, 903, L31

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

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

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

Chan M. H., Del Popolo A., 2020, MNRAS, 492, 5865

Cover T. M., Thomas J. A., 1991, Elements of Information Theory, 2nd edn. Wiley, Hoboken, New Jersey
not to say either that data could not plausibly have been generated by MOND or that it could plausibly have been generated under another hypothesis, as only MOND currently has sufficient predictivity for a test of this precision. We look to future SR studies with more data to establish the functional form of the RAR – if it exists – definitively.

## 6 CONCLUSION

The RAR has become central to debates about the mass discrepancy problem on astrophysical scales. Its tightness and regularity have been used to argue for a violation of Newtonian gravity in accordance with MOND, but the functions used to fit the data have been constructed to conform to this theory. As the first detailed application of the brand-new technique of ESR, we rank objectively *all* simple functions in terms of their aptitude for describing the SPARC RAR. We employ the minimum description length principle to trade accuracy with simplicity and hence perform model selection, and calibrate our method on mock MOND data generated both with and without the external field effect (EFE). Our conclusions are as follows:

(i) ESR discovers functions which are better descriptions, in both accuracy and simplicity and for both observed and simulated data, than MOND functions or a double power law.

(ii) While the majority of best-fitting functions on the SPARC data recover $g_{\text{obs}} \propto g_{\text{bar}}$ at high accelerations, not all have a best-fit coefficient of proportionality near unity. Thus the Newtonian limit is not clearly evidenced.

(iii) The SPARC data do not prefer functions with the deep-MOND limit of $g_{\text{obs}} \propto \sqrt{g_{\text{bar}}}$ as $g_{\text{bar}} \to 0$. Instead, we find that functions with $g_{\text{bar}} \to \text{const}$ typically compress the data more efficiently, albeit with considerable uncertainty.

(iv) SPARC-like mock data generated assuming the MONDian RAR interpolating function do not unambiguously recover that function. Moreover, many of the best functions for those mock data have $g_{\text{obs}} \approx 0.64 g_{\text{bar}}$ rather than $g_{\text{obs}} = g_{\text{bar}}$ at high $g_{\text{bar}}$, and most do not have a deep-MOND limit at all.

(v) The EFE in AQUAL greatly increases the logarithmic slope of the best-fitting functions at the low-$g_{\text{bar}}$ end of the data, but does not appreciably impact the limiting slope at $g_{\text{bar}} \to 0$. Incorporating the EFE in the mock data produces more generally similar results to the real data, so our analysis (within the MOND paradigm) hints at it.

(vi) We conclude that the data have too small a dynamic range (and too large uncertainties) to unambiguously favour MOND even if it is in fact generating the data. The SPARC RAR alone therefore does not supporting that theory unambiguously. The best prospect for improving this situation is to increase the acceleration range of the data, e.g. using stacked weak lensing at low $g_{\text{bar}}$ and groups and clusters at high $g_{\text{bar}}$.

(vii) Our results are a function of the maximum complexity of equation considered. Future symbolic regression algorithms – exhaustive or non-exhaustive – will reach the true description length minimum, and hence uncover the optimal functional representation of the RAR and determine whether the relation implies novel law-like gravitational behaviour.

ESR provides for the first time a guaranteed complete search through functional parameter space, making it the ideal tool to determine the analytic form of observed relations, extract physics from data theory-agnostically, and create fitting functions. We make the ESR and RAR codes, full function sets and the best 50 functions for each data set we consider publicly available to facilitate future applications.

## ACKNOWLEDGEMENTS

We thank Kyu-Hyun Chae, Andrei Constantin, Miles Cranmer, Mario Figueiredo, Gianluca Gregori, Thomas Harvey, Mark Kotanchek, Federico Lelli, Stacy McGaugh, Andre Lukas, Richard Stiskalek, and Tariq Yasin for useful inputs and discussion.

HD is supported by a Royal Society University Research Fellowship (grant no. 211046). DJB is supported by the Simons Collaboration on ‘Learning the Universe’ and was supported by STFC and Oriel College Oxford. PGF acknowledges support from European Research Council Grant No: 693024 and the Beecroft Trust.

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

This work used the DiRAC Complexity and DiRAC@Durham facilities, operated by the University of Leicester IT Services and Institute for Computational Cosmology, which form part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment is funded by Business Innovation & Skills (BIS) National E-Infrastructure capital grants ST/K000373/1, ST/P002293/1, ST/R002371/1, and ST/S002502/1, STFC DiRAC Operations grant ST/K0003259/1, and Durham University and Science and Technology Facilities Council (STFC) operations grant ST/R000832/1. DiRAC is part of the National E-Infrastructure.

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

## 7 DATA AVAILABILITY

The code and data associated with ESR and its application to the RAR are released at https://github.com/DeaglanBartlett/ESR and in Bartlett, Desmond & Ferreira (2022a). The SPARC data is available at http://astroweb.cwru.edu/SPARC. Other data may be shared on request to the corresponding authors.

## REFERENCES

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

Bartlett D. J., Desmond H., Ferreira P. G., 2022a, Exhaustive Symbolic Regression Function Sets

Bartlett D. J., Desmond H., Ferreira P. G., 2022b, preprint (arXiv:2211.11461)

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

Berezhiani L., Khoury J., 2015, Phys. Rev. D, 92, 103510

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

Blanchet L., Le Tiec A., 2008, Phys. Rev. D, 78, 024031

Brouwer M. M. et al., 2021, A&A, 650, A113

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

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

Chae K.-H., Bernardi M., Sheth R. K., Gong I.-T., 2019, ApJ, 877, 18

Chae K.-H., Bernardi M., Domínguez Sánchez H., Sheth R. K., 2020a, ApJ, 903, L31

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

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

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

Chan M. H., Del Popolo A., 2020, MNRAS, 492, 5865

Cover T. M., Thomas J. A., 1991, Elements of Information Theory, 2nd edn. Wiley, Hoboken, New Jersey
Cranmer M., Sanchez-Gonzalez A., Battaglia P., Xu R., Cranmer K., Spergel D., Ho S., 2020, in Larochelle H., Ranzato M., Hadsell R., Balcan M. F., Lin H., eds, *Advances in Neural Information Processing Systems*, vol. 33. Curran Associates, Inc., p. 17429

Desmond H., 2017, MNRAS, 464, 4160

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

Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161

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

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

Famaey B., Khoury J., Penco R., 2018, J. Cosmology Astropart. Phys., 2018, 038

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

Gopika K., Desai S., 2021, Phys. Dark Universe, 33, 100874

Grudić M. Y., Boylan-Kolchin M., Faucher-Giguère C.-A., Hopkins P. F., 2020, MNRAS, 496, L127

Grunwald P., 2007, The Minimum Description Length Principle. MIT Press, Massachusetts

Grünwald P., Roos T., 2019, *International Journal of Mathematics for Industry*, 11, 1930001

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

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

Kammerer L., Kronberger G., Burlacu B., Winkler S. M., Kommenda M., Affenzeller M., 2020, *Genetic Programming Theory and Practice XVII*. Springer, Cham, p. 79

Keller B. W., Wadsley J. W., 2017, ApJ, 835, L17

La Cava W., Orzechowski P., Burlacu B., Olivetti de França F., Virgolin M., Jin Y., Kommenda M., Moore J. H., 2021, in Vanschoren J., Yeung S., eds, *Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks*, vol. 1

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

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

Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3

Ludlow A. D. et al., 2017, Phys. Rev. Lett., 118, 161103

Macciò A. V., Stinson G., Brook C. B., Wadsley J., Couchman H. M. P., Shen S., Gibson B. K., Quinn T., 2012, ApJ, 744, L9

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

McGaugh S. S., 2004, ApJ, 609, 652

McGaugh S. S., Wolf J., 2010, ApJ, 722, 248

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

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

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

Milgrom M., 2010, MNRAS, 403, 886

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

Navarro J. F., Benítez-Llambay A., Fattahi A., Frenk C. S., Ludlow A. D., Oman K. A., Schaller M., Theuns T., 2017, MNRAS, 471, 1841

Oman K. A., Brouwer M. M., Ludlow A. D., Navarro J. F., 2020, preprint (arXiv:2006.06700)

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

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

Pradyumna S., Desai S., 2021, Phys. Dark Universe, 33, 100854

Rissanen J., 1978, Automatica, 14, 465

Rivero D., Fernandez-Blanco E., 2022, *Expert Systems with Applications*, 198, 116712

Sanders R. H., 1990, A&A Rev., 2, 1

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

Tenneti A., Mao Y.-Y., Croft R. A. C., Di Matteo T., Kosowsky A., Zago F., Zentner A. R., 2018, MNRAS, 474, 3125

Tian Y., Umetsu K., Ko C.-M., Donahue M., Chiu I. N., 2020, ApJ, 896, 70

Worm T., Chiu K., 2013, in Proceedings of the 15th Annual Conference on Genetic and Evolutionary Computation. GECCO '13. Association for Computing Machinery, New York, p. 1021

Zonoozi A. H., Lieberz P., Banik I., Haghi H., Kroupa P., 2021, MNRAS, 506, 5468

## APPENDIX A: A NOTE ON LIKELIHOODS

### A1 Correlation of measurements and their uncertainties

The likelihood of Section 3.2 treats the uncertainties induced by $\Upsilon$, $D$, and $i$ as statistical and uncorrelated between points. This is clearly incorrect: a scattering up of $\Upsilon$, for example, causes a coherent increase in $g_{bar}$ across the rotation curve, generating off-diagonal elements in the covariance matrix. A better approach is therefore to calculate the covariance matrix via Monte Carlo. For a given galaxy, one would independently sample $\Upsilon_{gas}$, $\Upsilon_{disc}$, $\Upsilon_{bulge}$, $D$, and $i$ many times from their assumed-Gaussian prior distributions. One would then generate the corresponding $g_{obs}$ and $g_{bar}$ values (now vectors across the rotation curve of a given galaxy), scatter them by their statistical uncertainty ($g_{bar}$ has only the small $\delta L_{3,6}$ term), and calculate the covariance matrices $\Sigma_{obs}$ and $\Sigma_{bar}$.

Let us define $u \equiv \log(g_{bar})$ and $v \equiv \log(g_{obs})$ such that $u$ and $v$ are vectors containing $u$ and $v$ for all galaxies at all measured points along the rotation curve. We denote the 'true' values with a superscript t and the observed values without a superscript. We assume that $u^t$ and $v^t$ are drawn from multivariate Gaussian distributions with covariance matrices $\Sigma_u$ and $\Sigma_v$, respectively. Since $u^t$ and $v^t$ are statistically independent random variables for fixed $u$ and $v$, the likelihood of a given set of points in the $(u^t, v^t)$ plane would then be

$$ \mathcal{L}(u, v) = \frac{1}{\sqrt{|2\pi\Sigma_u|} \sqrt{|2\pi\Sigma_v|}} \times \exp\left(-\frac{1}{2}(u^t - u)^T \Sigma_u^{-1} (u^t - u)\right) \times \exp\left(-\frac{1}{2}(v^t - v)^T \Sigma_v^{-1} (v^t - v)\right). \quad (A1) $$

The equivalent of equation (2) for a vector-valued function $f$, $v = f(u)$, is then obtained by marginalising over $u^t$,

$$ \mathcal{L}(v) = \int \mathcal{L}(u^t, f(u^t)) du^t. \quad (A2) $$

Taylor expanding $f(u^t)$ about $u$,

$$ f(u^t) \approx f(u) + D(u^t - u), \quad (A3) $$

where $D_{ij} \equiv \partial_j f_i|_u$, yields

$$ \mathcal{L}(v) = \frac{1}{\sqrt{|2\pi\Sigma_{tot}|}} \times \exp\left(-\frac{1}{2}(v - f(u))^T \Sigma_{tot}^{-1} (v - f(u))\right), \quad (A4) $$

where $\Sigma_{tot}$ is the total covariance matrix defined by

$$ \Sigma_{tot} \equiv \Sigma_v + D \Sigma_u D^T. \quad (A5) $$

This is equivalent to equation (3) if $u$ and $v$ have only a single element.

In principle, the assumption that the joint probability distributions of $u^t$ and $v^t$ are a multivariate Gaussian is unnecessary as the full empirical distributions are generated by the Monte Carlo sampling described above. Using this directly would enable a loss function which, while still independent between galaxies, fully encapsulates the correlated, non-Gaussian structure of each galaxy's measurements, and hence provides a more accurate description of the expected probability distributions of $g_{bar}$ and $g_{obs}$. This model would also allow for more accurate mock data generation as discussed in Section 5.
Cranmer M., Sanchez-Gonzalez A., Battaglia P., Xu R., Cranmer K., Spergel D., Ho S., 2020, in Larochelle H., Ranzato M., Hadsell R., Balcan M. F., Lin H., eds, *Advances in Neural Information Processing Systems*, vol. 33. Curran Associates, Inc., p. 17429

Desmond H., 2017, MNRAS, 464, 4160

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

Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., Battye R. A., Booth C. M., 2010, MNRAS, 405, 2161

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

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

Famaey B., Khoury J., Penco R., 2018, J. Cosmology Astropart. Phys., 2018, 038

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

Gopika K., Desai S., 2021, Phys. Dark Universe, 33, 100874

Grudić M. Y., Boylan-Kolchin M., Faucher-Giguère C.-A., Hopkins P. F., 2020, MNRAS, 496, L127

Grunwald P., 2007, The Minimum Description Length Principle. MIT Press, Massachusetts

Grünwald P., Roos T., 2019, *International Journal of Mathematics for Industry*, 11, 1930001

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

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

Kammerer L., Kronberger G., Burlacu B., Winkler S. M., Kommenda M., Affenzeller M., 2020, *Genetic Programming Theory and Practice XVII*. Springer, Cham, p. 79

Keller B. W., Wadsley J. W., 2017, ApJ, 835, L17

La Cava W., Orzechowski P., Burlacu B., Olivetti de França F., Virgolin M., Jin Y., Kommenda M., Moore J. H., 2021, in Vanschoren J., Yeung S., eds, *Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks*, vol. 1

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

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

Li P., Lelli F., McGaugh S., Schombert J., 2018, A&A, 615, A3

Ludlow A. D. et al., 2017, Phys. Rev. Lett., 118, 161103

Macciò A. V., Stinson G., Brook C. B., Wadsley J., Couchman H. M. P., Shen S., Gibson B. K., Quinn T., 2012, ApJ, 744, L9

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

McGaugh S. S., 2004, ApJ, 609, 652

McGaugh S. S., Wolf J., 2010, ApJ, 722, 248

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

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

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

Milgrom M., 2010, MNRAS, 403, 886

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

Navarro J. F., Benítez-Llambay A., Fattahi A., Frenk C. S., Ludlow A. D., Oman K. A., Schaller M., Theuns T., 2017, MNRAS, 471, 1841

Oman K. A., Brouwer M. M., Ludlow A. D., Navarro J. F., 2020, preprint (arXiv:2006.06700)

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

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

Pradyumna S., Desai S., 2021, Phys. Dark Universe, 33, 100854

Rissanen J., 1978, Automatica, 14, 465

Rivero D., Fernandez-Blanco E., 2022, *Expert Systems with Applications*, 198, 116712

Sanders R. H., 1990, A&A Rev., 2, 1

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

Tenneti A., Mao Y.-Y., Croft R. A. C., Di Matteo T., Kosowsky A., Zago F., Zentner A. R., 2018, MNRAS, 474, 3125

Tian Y., Umetsu K., Ko C.-M., Donahue M., Chiu I. N., 2020, ApJ, 896, 70

Worm T., Chiu K., 2013, in Proceedings of the 15th Annual Conference on Genetic and Evolutionary Computation. GECCO '13. Association for Computing Machinery, New York, p. 1021

Zonoozi A. H., Lieberz P., Banik I., Haghi H., Kroupa P., 2021, MNRAS, 506, 5468

## APPENDIX A: A NOTE ON LIKELIHOODS

### A1 Correlation of measurements and their uncertainties

The likelihood of Section 3.2 treats the uncertainties induced by $\Upsilon$, $D$, and $i$ as statistical and uncorrelated between points. This is clearly incorrect: a scattering up of $\Upsilon$, for example, causes a coherent increase in $g_{bar}$ across the rotation curve, generating off-diagonal elements in the covariance matrix. A better approach is therefore to calculate the covariance matrix via Monte Carlo. For a given galaxy, one would independently sample $\Upsilon_{gas}$, $\Upsilon_{disc}$, $\Upsilon_{bulge}$, $D$, and $i$ many times from their assumed-Gaussian prior distributions. One would then generate the corresponding $g_{obs}$ and $g_{bar}$ values (now vectors across the rotation curve of a given galaxy), scatter them by their statistical uncertainty ($g_{bar}$ has only the small $\delta L_{3,6}$ term), and calculate the covariance matrices $\Sigma_{obs}$ and $\Sigma_{bar}$.

Let us define $u \equiv \log(g_{bar})$ and $v \equiv \log(g_{obs})$ such that $u$ and $v$ are vectors containing $u$ and $v$ for all galaxies at all measured points along the rotation curve. We denote the 'true' values with a superscript t and the observed values without a superscript. We assume that $u^t$ and $v^t$ are drawn from multivariate Gaussian distributions with covariance matrices $\Sigma_u$ and $\Sigma_v$, respectively. Since $u^t$ and $v^t$ are statistically independent random variables for fixed $u$ and $v$, the likelihood of a given set of points in the $(u^t, v^t)$ plane would then be

$$ \mathcal{L}(u, v) = \frac{1}{\sqrt{|2\pi\Sigma_u|}} \frac{1}{\sqrt{|2\pi\Sigma_v|}} \times \exp\left(-\frac{1}{2}(u^t - u)^T \Sigma_u^{-1} (u^t - u)\right) \times \exp\left(-\frac{1}{2}(v^t - v)^T \Sigma_v^{-1} (v^t - v)\right). \quad (A1) $$

The equivalent of equation (2) for a vector-valued function $f$, $v = f(u)$, is then obtained by marginalising over $u^t$,

$$ \mathcal{L}(v) = \int \mathcal{L}(u^t, f(u^t)) du^t. \quad (A2) $$

Taylor expanding $f(u^t)$ about $u$,

$$ f(u^t) \approx f(u) + D(u^t - u), \quad (A3) $$

where $D_{ij} \equiv \partial_j f_i|_u$, yields

$$ \mathcal{L}(v) = \frac{1}{\sqrt{|2\pi\Sigma_{tot}|}} \times \exp\left(-\frac{1}{2}(v - f(u))^T \Sigma_{tot}^{-1} (v - f(u))\right), \quad (A4) $$

where $\Sigma_{tot}$ is the total covariance matrix defined by

$$ \Sigma_{tot} \equiv \Sigma_v + D \Sigma_u D^T. \quad (A5) $$

This is equivalent to equation (3) if $u$ and $v$ have only a single element.

In principle, the assumption that the joint probability distributions of $u^t$ and $v^t$ are a multivariate Gaussian is unnecessary as the full empirical distributions are generated by the Monte Carlo sampling described above. Using this directly would enable a loss function which, while still independent between galaxies, fully encapsulates the correlated, non-Gaussian structure of each galaxy's measurements, and hence provides a more accurate description of the expected probability distributions of $g_{bar}$ and $g_{obs}$. This model would also allow for more accurate mock data generation as discussed in Section 5.
## A2 Sampling nuisance parameters from their posteriors

All the methods described above sample $\Upsilon, D$, and $i$, and the true $g_{bar}$ values $u^t$, from their priors, i.e. without adjusting them to maximize agreement with the function being fitted. In principle, the better procedure is to constrain these nuisance parameters jointly with any parameters of the function, assuming that function in the fit. Even the relatively simple case where $u^t$ is drawn from the prior but $\Upsilon, D$, and $i$ from the posterior is however impractical as it requires optimization in a parameter space of dimension $147 \times 3 + p$ if the three $\Upsilon$s are coupled (or only one is varied), and $147 \times 5 + p$ if they are varied separately, where $p$ is the number of free parameters in the function. This would naturally account for correlations between the $g_{bar}$ and $g_{obs}$ measurements induced by variations in $\Upsilon, D$, and $i$, and thus provides a more accurate but expensive alternative to the model of Section A1. This approach will be applied for the first time to the RAR, outside the context of SR, in upcoming work (Desmond 2023, in prep).

The situation is complicated further by inference of $u^t$. At the top level of the hierarchical model for predicting $v$ from $u$ are both the parameters of the model, $\theta$, and $u^t$. If one wanted to find the true maximum likelihood point, one should maximize equation (A1) for $u^t$ as well as $\theta$, instead of equation (A2) for $\theta$ alone.

To do this, we start by noting that maximising the likelihood in equation (A1) is equivalent to minimising

$$h(u^t, v) = \frac{1}{2}(f(u^t) - v)^T \Sigma_v^{-1} (f(u^t) - v) + \frac{1}{2}(u^t - u)^T \Sigma_u^{-1} (u^t - u). \quad (A6)$$

Using equation (A3) and minimising $h$ with respect to $u^t$, we find the maximum-likelihood point, $\hat{u}^t$, to be

$$\hat{u}^t = u + \sum_u D^T \Sigma_{tot}^{-1} (v - f(u)), \quad (A7)$$

and we must therefore maximize

$$\mathcal{L}(\hat{u}^t, v) = \frac{1}{\sqrt{|2\pi\Sigma_u|}} \frac{1}{\sqrt{|2\pi\Sigma_v|}} \times \exp\left(-\frac{1}{2}(v - f(u))^T \Sigma_{tot}^{-1} (v - f(u))\right), \quad (A8)$$

where $\Sigma_{tot}^{-1}$ is defined in equation (A5). This will not yield the same result as maximising equation (A4) since equation (A8) does not contain the normalization term which penalizes large gradients. When applied to mock data generated assuming $v^t$ is linearly related

to $u^t$, we find that equation (A4) (the commonly used expression, and the one we adopt in Section 3.2) induces a small bias that equation (A8) does not. This is the reason why the best-fit RAR and generalized RAR IFs do not precisely match the generating function in Table 2, or the Simple IF + EFE the generating function in Table 3. We should expect all best-fitting parameter values to be likewise slightly biased by the use of equation (2).

In practice, it is challenging to maximize equation (A1) since the lack of the gradient-penalizing determinant means that naively using equation (A8) can prefer functions with diverging gradients at at least one point in the domain of $u$. This breaks the linearity assumption in equation (A3), making the result untrustworthy. Instead, one should numerically solve the full optimization problem without Taylor-expanding $f$. This involves solving a root-finding problem for each trial $\theta$ during the optimization of the function's parameters. Although feasible for a handful of functions, this is computationally impractical for the full set of ESR functions. We therefore use the simpler likelihood here, calibrating our results using mock data to measure the magnitude of the bias. We defer further discussion of this important issue to future work.

## APPENDIX B: LIMITING SLOPES OF ESR FUNCTIONS

Here we provide the analytic low-$g_{bar}$ ($s_-$) and high-$g_{bar}$ ($s_+$) logarithmic slopes of the top 10 functions generated using ESR up to complexity 9. We give the slopes for the observed SPARC data in Table B1, the mock data assuming the RAR IF in Table B2, and the mock data assuming the Simple IF with EFE in Table B3.
## A2 Sampling nuisance parameters from their posteriors

All the methods described above sample $\Upsilon, D$, and $i$, and the true $g_{bar}$ values $u^t$, from their priors, i.e. without adjusting them to maximize agreement with the function being fitted. In principle, the better procedure is to constrain these nuisance parameters jointly with any parameters of the function, assuming that function in the fit. Even the relatively simple case where $u^t$ is drawn from the prior but $\Upsilon, D$, and $i$ from the posterior is however impractical as it requires optimization in a parameter space of dimension $147 \times 3 + p$ if the three $\Upsilon$s are coupled (or only one is varied), and $147 \times 5 + p$ if they are varied separately, where $p$ is the number of free parameters in the function. This would naturally account for correlations between the $g_{bar}$ and $g_{obs}$ measurements induced by variations in $\Upsilon, D$, and $i$, and thus provides a more accurate but expensive alternative to the model of Section A1. This approach will be applied for the first time to the RAR, outside the context of SR, in upcoming work (Desmond 2023, in prep).

The situation is complicated further by inference of $u^t$. At the top level of the hierarchical model for predicting $v$ from $u$ are both the parameters of the model, $\theta$, and $u^t$. If one wanted to find the true maximum likelihood point, one should maximize equation (A1) for $u^t$ as well as $\theta$, instead of equation (A2) for $\theta$ alone.

To do this, we start by noting that maximising the likelihood in equation (A1) is equivalent to minimising

$$h(u^t, v) = \frac{1}{2}(f(u^t) - v)^T \Sigma_v^{-1} (f(u^t) - v) + \frac{1}{2}(u^t - u)^T \Sigma_u^{-1} (u^t - u). \quad (A6)$$

Using equation (A3) and minimising $h$ with respect to $u^t$, we find the maximum-likelihood point, $\hat{u}^t$, to be

$$\hat{u}^t = u + \sum_u D^T \Sigma_{\text{tot}}^{-1} (v - f(u)), \quad (A7)$$

and we must therefore maximize

$$\mathcal{L}(\hat{u}^t, v) = \frac{1}{\sqrt{|2\pi\Sigma_u|}} \frac{1}{\sqrt{|2\pi\Sigma_v|}} \times \exp\left(-\frac{1}{2}(v - f(u))^T \Sigma_{\text{tot}}^{-1} (v - f(u))\right), \quad (A8)$$

where $\Sigma_{\text{tot}}^{-1}$ is defined in equation (A5). This will not yield the same result as maximising equation (A4) since equation (A8) does not contain the normalization term which penalizes large gradients. When applied to mock data generated assuming $v^t$ is linearly related

to $u^t$, we find that equation (A4) (the commonly used expression, and the one we adopt in Section 3.2) induces a small bias that equation (A8) does not. This is the reason why the best-fit RAR and generalized RAR IFs do not precisely match the generating function in Table 2, or the Simple IF + EFE the generating function in Table 3. We should expect all best-fitting parameter values to be likewise slightly biased by the use of equation (2).

In practice, it is challenging to maximize equation (A1) since the lack of the gradient-penalizing determinant means that naively using equation (A8) can prefer functions with diverging gradients at at least one point in the domain of $u$. This breaks the linearity assumption in equation (A3), making the result untrustworthy. Instead, one should numerically solve the full optimization problem without Taylor-expanding $f$. This involves solving a root-finding problem for each trial $\theta$ during the optimization of the function's parameters. Although feasible for a handful of functions, this is computationally impractical for the full set of ESR functions. We therefore use the simpler likelihood here, calibrating our results using mock data to measure the magnitude of the bias. We defer further discussion of this important issue to future work.

## APPENDIX B: LIMITING SLOPES OF ESR FUNCTIONS

Here we provide the analytic low-$g_{bar}$ ($s_-$) and high-$g_{bar}$ ($s_+$) logarithmic slopes of the top 10 functions generated using ESR up to complexity 9. We give the slopes for the observed SPARC data in Table B1, the mock data assuming the RAR IF in Table B2, and the mock data assuming the Simple IF with EFE in Table B3.
<table><caption>**Table B1.** Functional forms and limiting slopes of the ten best functions found by ESR applied to the SPARC data. The functions in the second column give the fitted $y = g_{\text{obs}}/10^{-10}$ m s$^{-2}$ for input $x = g_{\text{bar}}/10^{-10}$ m s$^{-2}$. The low-acceleration slope is $\lim_{x \to 0^+} \text{d log y/d log x}$ (denoted $s_-$ in the text), and the high-acceleration slope is similarly defined but for $x \to \infty$ ($s_+$). {$\theta_i$} are real parameters fitted to the data to maximize the likelihood (see Table 1). Slopes given without conditions are valid $\forall \theta_i$. For comparison, the MOND prediction without the external field effect is $s_- = 1/2$ and $s_+ = 1$.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th colspan="2">Low-acceleration slope</th><th colspan="2">High-acceleration slope</th></tr><tr><th>Value</th><th>Condition</th><th>Value</th><th>Condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> (|θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup> + x)</td><td>0</td><td></td><td>θ<sub>2</sub></td><td>θ<sub>2</sub> ≥ 1</td></tr><tr><td>2</td><td>|θ<sub>1</sub>|<sup>x</sup> + θ<sub>0</sub>|<sup>θ<sub>2</sub></sup> + x</td><td>0</td><td></td><td>1</td><td>otherwise</td></tr><tr><td>3</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub> - x|<sup>θ<sub>2</sub> - θ<sub>3</sub></sup></sub></td><td>0</td><td></td><td>1</td><td>(|θ<sub>1</sub>| ≤ 1) or (|θ<sub>1</sub>| > 1 and θ<sub>2</sub> ≤ 0)</td></tr><tr><td>4</td><td>|θ<sub>0</sub>(θ<sub>1</sub> + x)|<sup>θ<sub>2</sub></sup> + x</td><td>0</td><td></td><td>∞</td><td>otherwise</td></tr><tr><td>5</td><td>|θ<sub>0</sub> - |θ<sub>1</sub> - x|<sup>θ<sub>2</sub></sup>|<sup>θ<sub>3</sub></sup></td><td>0</td><td></td><td>∞</td><td>|θ<sub>0</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>6</td><td>√x exp (<sup>|θ<sub>0</sub>+x|<sup>θ<sub>1</sub></sup></sup> / 2)</td><td>1/2</td><td></td><td>-∞</td><td>0 < |θ<sub>0</sub>| < 1 and θ<sub>2</sub> > 0</td></tr><tr><td>7</td><td>(<sup>|θ<sub>0</sub>|<sup>x</sup></sup> / x)<sup>θ<sub>1</sub></sup> + x</td><td>-θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≥ -1<br>otherwise</td><td>1<br>-θ<sub>1</sub><br>1 + 0<sup>θ<sub>1</sub></sup> θ<sub>1</sub>(-∞)</td><td>(|θ<sub>0</sub>| > 1 and θ<sub>1</sub> ≤ 0) or (0 < |θ<sub>0</sub>| < 1 and θ<sub>1</sub> ≥ 0) or (|θ<sub>0</sub>| = 1 and θ<sub>1</sub> ≥ -1)<br>|θ<sub>0</sub>| = 1 and θ<sub>1</sub> < -1<br>θ<sub>0</sub> = 0<br>otherwise</td></tr><tr><td>8</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>0</td><td></td><td>1</td><td>-</td></tr><tr><td>9</td><td>|θ<sub>0</sub> + 1/√x|<sup>θ<sub>1</sub></sup></td><td>-θ<sub>1</sub>/4</td><td></td><td>0</td><td>-</td></tr><tr><td>10</td><td>(√x + 1/√x)<sup>θ<sub>0</sub></sup> + x</td><td>-θ<sub>0</sub><br>1</td><td>θ<sub>0</sub> ≥ -1<br>otherwise</td><td>θ<sub>0</sub>/2<br>1</td><td>θ<sub>0</sub> ≥ 2<br>otherwise</td></tr></tbody></table>

<table><caption>**Table B2.** As Table B1 but for the mock data generated from the RAR IF.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th colspan="2">Low-acceleration slope</th><th colspan="2">High-acceleration slope</th></tr><tr><th>Value</th><th>Condition</th><th>Value</th><th>Condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> + θ<sub>1</sub>x + √x</td><td>0</td><td>-</td><td>1</td><td>-</td></tr><tr><td>2</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>0</td><td>-</td><td>1</td><td>-</td></tr><tr><td>3</td><td>θ<sub>0</sub>x + x<sup>θ<sub>1</sub></sup></td><td>θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≤ 1<br>otherwise</td><td>θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≥ 1<br>otherwise</td></tr><tr><td>4</td><td>√x exp (<sup>xθ<sub>0</sub></sup> / 2)</td><td>1/2<br>-∞</td><td>θ<sub>0</sub> ≥ 0<br>otherwise</td><td>1/2<br>∞</td><td>θ<sub>0</sub> ≤ 0<br>otherwise</td></tr><tr><td>5</td><td>(θ<sub>0</sub> + x) (θ<sub>1</sub> + 1/√x)</td><td>-1/2</td><td>-</td><td>1</td><td>-</td></tr><tr><td>6</td><td>1/√(|θ<sub>0</sub> + 1/√x|) + x</td><td>1/2</td><td>-</td><td>1</td><td>-</td></tr><tr><td>7</td><td>(x|θ<sub>0</sub>|)(x|θ<sub>1</sub>|)<sup>θ<sub>2</sub></sup></td><td>∞<br>1<br>0</td><td>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> < 0<br>θ<sub>2</sub> = 0<br>otherwise</td><td>∞<br>0<br>1</td><td>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> > 0<br>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> < 0<br>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> = 0</td></tr><tr><td>8</td><td>θ<sub>0</sub>x + |θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup></td><td>0</td><td>-<br>-<br>0</td><td>θ<sub>2</sub><br>1</td><td>θ<sub>2</sub> ≥ 1<br>otherwise</td></tr><tr><td>9</td><td>x (|θ<sub>0</sub> - x|<sup>θ<sub>1</sub> - θ<sub>2</sub></sup>)</td><td>1</td><td>-<br>-</td><td>θ<sub>1</sub> + 1<br>1</td><td>θ<sub>1</sub> ≥ 0<br>otherwise</td></tr><tr><td>10</td><td>(θ<sub>0</sub> - x) (θ<sub>1</sub> - x<sup>θ<sub>2</sub></sup>)</td><td>θ<sub>2</sub><br>0</td><td>θ<sub>2</sub> ≤ 0<br>otherwise</td><td>θ<sub>2</sub> + 1<br>1</td><td>θ<sub>2</sub> ≥ 0<br>otherwise</td></tr></tbody></table>
<table><caption>**Table B1.** Functional forms and limiting slopes of the ten best functions found by ESR applied to the SPARC data. The functions in the second column give the fitted $y = g_{\text{obs}}/10^{-10}$ m s$^{-2}$ for input $x = g_{\text{bar}}/10^{-10}$ m s$^{-2}$. The low-acceleration slope is $\lim_{x \to 0^+} \text{d log y/d log x}$ (denoted $s_-$ in the text), and the high-acceleration slope is similarly defined but for $x \to \infty$ ($s_+$). {$\theta_i$} are real parameters fitted to the data to maximize the likelihood (see Table 1). Slopes given without conditions are valid $\forall \theta_i$. For comparison, the MOND prediction without the external field effect is $s_- = 1/2$ and $s_+ = 1$.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th colspan="2">Low-acceleration slope</th><th colspan="2">High-acceleration slope</th></tr><tr><th>Value</th><th>Condition</th><th>Value</th><th>Condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> (|θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup> + x)</td><td>0</td><td></td><td>θ<sub>2</sub></td><td>θ<sub>2</sub> ≥ 1</td></tr><tr><td>2</td><td>|θ<sub>1</sub>|<sup>x</sup> + θ<sub>0</sub>|<sup>θ<sub>2</sub></sup> + x</td><td>0</td><td></td><td>1</td><td>otherwise</td></tr><tr><td>3</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub> - x|<sup>θ<sub>2</sub> - θ<sub>3</sub></sup></sup></td><td>0</td><td></td><td>1</td><td>(|θ<sub>1</sub>| ≤ 1) or (|θ<sub>1</sub>| > 1 and θ<sub>2</sub> ≤ 0)</td></tr><tr><td>4</td><td>|θ<sub>0</sub>(θ<sub>1</sub> + x)|<sup>θ<sub>2</sub></sup> + x</td><td>0</td><td></td><td>∞</td><td>otherwise</td></tr><tr><td>5</td><td>|θ<sub>0</sub> - |θ<sub>1</sub> - x|<sup>θ<sub>2</sub></sup>|<sup>θ<sub>3</sub></sup></td><td>0</td><td></td><td>∞</td><td>|θ<sub>0</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>6</td><td>√x exp (<sup>|θ<sub>0</sub>+x|<sup>θ<sub>1</sub></sup></sup> / 2)</td><td>1/2</td><td></td><td>-∞</td><td>0 < |θ<sub>0</sub>| < 1 and θ<sub>2</sub> > 0</td></tr><tr><td>7</td><td>(<sup>|θ<sub>0</sub>|<sup>x</sup></sup> / x)<sup>θ<sub>1</sub></sup> + x</td><td>-θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≥ -1<br>otherwise</td><td>1<br>-θ<sub>1</sub><br>1 + 0<sup>θ<sub>1</sub></sup> θ<sub>1</sub>(-∞)</td><td>(|θ<sub>0</sub>| > 1 and θ<sub>1</sub> ≤ 0) or (0 < |θ<sub>0</sub>| < 1 and θ<sub>1</sub> ≥ 0) or (|θ<sub>0</sub>| = 1 and θ<sub>1</sub> ≥ -1)<br>|θ<sub>0</sub>| = 1 and θ<sub>1</sub> < -1<br>θ<sub>0</sub> = 0<br>otherwise</td></tr><tr><td>8</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>0</td><td></td><td>1</td><td>-</td></tr><tr><td>9</td><td>|θ<sub>0</sub> + 1/√x|<sup>θ<sub>1</sub></sup></td><td>-θ<sub>1</sub>/4</td><td></td><td>0</td><td>-</td></tr><tr><td>10</td><td>(√x + 1/√x)<sup>θ<sub>0</sub></sup> + x</td><td>-θ<sub>0</sub><br>1</td><td>θ<sub>0</sub> ≥ -1<br>otherwise</td><td>θ<sub>0</sub>/2<br>1</td><td>θ<sub>0</sub> ≥ 2<br>otherwise</td></tr></tbody></table>

<table><caption>**Table B2.** As Table B1 but for the mock data generated from the RAR IF.</caption><thead><tr><th rowspan="2">Rank</th><th rowspan="2">Function</th><th colspan="2">Low-acceleration slope</th><th colspan="2">High-acceleration slope</th></tr><tr><th>Value</th><th>Condition</th><th>Value</th><th>Condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> + θ<sub>1</sub>x + √x</td><td>0</td><td>-</td><td>1</td><td>-</td></tr><tr><td>2</td><td>√(|θ<sub>0</sub> + x| + θ<sub>1</sub>x)</td><td>0</td><td>-</td><td>1</td><td>-</td></tr><tr><td>3</td><td>θ<sub>0</sub>x + x<sup>θ<sub>1</sub></sup></td><td>θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≤ 1<br>otherwise</td><td>θ<sub>1</sub><br>1</td><td>θ<sub>1</sub> ≥ 1<br>otherwise</td></tr><tr><td>4</td><td>√x exp (<sup>xθ<sub>0</sub></sup> / 2)</td><td>1/2<br>-∞</td><td>θ<sub>0</sub> ≥ 0<br>otherwise</td><td>1/2<br>∞</td><td>θ<sub>0</sub> ≤ 0<br>otherwise</td></tr><tr><td>5</td><td>(θ<sub>0</sub> + x)(θ<sub>1</sub> + 1/√x)</td><td>-1/2</td><td>-</td><td>1</td><td>-</td></tr><tr><td>6</td><td>1/√(|θ<sub>0</sub> + 1/√x|) + x</td><td>1/2</td><td>-</td><td>1</td><td>-</td></tr><tr><td>7</td><td>(x|θ<sub>0</sub>|)(x|θ<sub>1</sub>|)<sup>θ<sub>2</sub></sup></td><td>∞<br>1<br>0</td><td>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> < 0<br>θ<sub>2</sub> = 0<br>otherwise</td><td>∞<br>0<br>1</td><td>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> > 0<br>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> < 0<br>|θ<sub>1</sub>| > 0 and θ<sub>2</sub> = 0</td></tr><tr><td>8</td><td>θ<sub>0</sub>x + |θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup></td><td>0</td><td>-<br>-<br>-</td><td>θ<sub>2</sub><br>1</td><td>θ<sub>2</sub> ≥ 1<br>otherwise</td></tr><tr><td>9</td><td>x (|θ<sub>0</sub> - x|<sup>θ<sub>1</sub> - θ<sub>2</sub></sup>)</td><td>1</td><td>-<br>-</td><td>θ<sub>1</sub> + 1<br>1</td><td>θ<sub>1</sub> ≥ 0<br>otherwise</td></tr><tr><td>10</td><td>(θ<sub>0</sub> - x)(θ<sub>1</sub> - x<sup>θ<sub>2</sub></sup>)</td><td>θ<sub>2</sub><br>0</td><td>θ<sub>2</sub> ≤ 0<br>otherwise</td><td>θ<sub>2</sub> + 1<br>1</td><td>θ<sub>2</sub> ≥ 0<br>otherwise</td></tr></tbody></table>
<table><caption>**Table B3.** As Table B1 but for the mock data generated using the Simple IF + EFE.</caption><thead><tr><th>Rank</th><th>Function</th><th>Value</th><th>Low-acceleration slope condition</th><th>Value</th><th>High-acceleration slope condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> + √(x<sup>2</sup> + 2x)</td><td>0</td><td>—</td><td>1</td><td>—</td></tr><tr><td>2</td><td>θ<sub>0</sub> + √(x|θ<sub>1</sub> + x|)</td><td>0</td><td>—</td><td>1</td><td>—</td></tr><tr><td>3</td><td>-|θ<sub>0</sub>|<sup>√x</sup> + θ<sub>1</sub> + x</td><td>0</td><td>—</td><td>1</td><td>|θ<sub>0</sub>| ≤ 1</td></tr><tr><td>4</td><td>(θ<sub>0</sub> - x)(θ<sub>1</sub> - x<sup>θ<sub>2</sub></sup>)</td><td>θ<sub>2</sub></td><td>θ<sub>2</sub> ≤ 0</td><td>1 + θ<sub>2</sub></td><td>otherwise</td></tr><tr><td>5</td><td>x<sup>θ<sub>0</sub></sup> - θ<sub>1</sub>(θ<sub>2</sub> - x)</td><td>θ<sub>0</sub></td><td>θ<sub>0</sub> ≤ 0</td><td>1</td><td>θ<sub>2</sub> ≥ 0</td></tr><tr><td>6</td><td>|θ<sub>0</sub> - x|<sup>θ<sub>1</sub></sup> - θ<sub>2</sub>x</td><td>0</td><td>otherwise</td><td>1</td><td>otherwise</td></tr><tr><td>7</td><td>x|θ<sub>0</sub>|<sup>-|θ<sub>1</sub>|<sup>x|θ<sub>2</sub></sup></sup></td><td>∞log|θ<sub>0</sub>|</td><td>|θ<sub>1</sub>| > 1 and θ<sub>2</sub> < 0</td><td>-∞</td><td>|θ<sub>0</sub>| > 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>8</td><td>x(θ<sub>0</sub> + |θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup>)</td><td>1</td><td>otherwise</td><td>1</td><td>0 < |θ<sub>0</sub>| < 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>9</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub>|<sup>x|θ<sub>2</sub></sup></sup> + x</td><td>-∞</td><td>|θ<sub>0</sub>| > 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> < 0</td><td>∞</td><td>θ<sub>0</sub> = 0 and |θ<sub>1</sub>| > 1</td></tr><tr><td>10</td><td>exp(θ<sub>0</sub> - 1/√x) + x</td><td>1</td><td>otherwise</td><td>1</td><td>otherwise</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.
<table><caption>**Table B3.** As Table B1 but for the mock data generated using the Simple IF + EFE.</caption><thead><tr><th>Rank</th><th>Function</th><th>Value</th><th>Low-acceleration slope condition</th><th>Value</th><th>High-acceleration slope condition</th></tr></thead><tbody><tr><td>1</td><td>θ<sub>0</sub> + √(x<sup>2</sup> + 2x)</td><td>0</td><td>–</td><td>1</td><td>–</td></tr><tr><td>2</td><td>θ<sub>0</sub> + √(x|θ<sub>1</sub> + x|)</td><td>0</td><td>–</td><td>1</td><td>–</td></tr><tr><td>3</td><td>−|θ<sub>0</sub>|<sup>√x</sup> + θ<sub>1</sub> + x</td><td>0</td><td>–</td><td>1</td><td>|θ<sub>0</sub>| ≤ 1</td></tr><tr><td>4</td><td>(θ<sub>0</sub> − x)(θ<sub>1</sub> − x<sup>θ<sub>2</sub></sup>)</td><td>θ<sub>2</sub></td><td>θ<sub>2</sub> ≤ 0</td><td>1 + θ<sub>2</sub></td><td>otherwise</td></tr><tr><td>5</td><td>x<sup>θ<sub>0</sub></sup> − θ<sub>1</sub>(θ<sub>2</sub> − x)</td><td>θ<sub>0</sub></td><td>θ<sub>0</sub> ≤ 0</td><td>1</td><td>θ<sub>2</sub> ≥ 0</td></tr><tr><td>6</td><td>|θ<sub>0</sub> − x|<sup>θ<sub>1</sub></sup> − θ<sub>2</sub>x</td><td>0</td><td>otherwise</td><td>1</td><td>otherwise</td></tr><tr><td>7</td><td>x|θ<sub>0</sub>|<sup>−|θ<sub>1</sub>|<sup>x</sup>θ<sub>2</sub></sup></td><td>∞log|θ<sub>0</sub>|</td><td>|θ<sub>1</sub>| > 1 and θ<sub>2</sub> < 0</td><td>−∞</td><td>|θ<sub>0</sub>| > 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>8</td><td>x(θ<sub>0</sub> + |θ<sub>1</sub> + x|<sup>θ<sub>2</sub></sup>)</td><td>1</td><td>1</td><td>∞</td><td>0 < |θ<sub>0</sub>| < 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> > 0</td></tr><tr><td>9</td><td>|θ<sub>0</sub>|<sup>|θ<sub>1</sub>|<sup>x</sup>θ<sub>2</sub></sup> + x</td><td>−∞</td><td>|θ<sub>0</sub>| > 1 and |θ<sub>1</sub>| > 1 and θ<sub>2</sub> < 0</td><td>∞</td><td>θ<sub>0</sub> = 0 and |θ<sub>1</sub>| > 1</td></tr><tr><td>10</td><td>exp(θ<sub>0</sub> − 1/√x) + x</td><td>1</td><td>1</td><td>1</td><td>otherwise</td></tr></tbody></table>

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