© ESO 2021

1. Introduction

The missing mass problem in astrophysical systems, from galaxies and galaxy clusters to the cosmic microwave background, is well established. Gravitational potentials are observed to be deeper than predicted from the visible matter distributions in Newtonian gravity. Early observations of this phenomenon include the velocity dispersion of galaxies in clusters (Zwicky 1933) and galactic rotation curves (Rubin & Ford 1970; Rubin et al. 1980; Bosma 1981). Both particle dark matter (DM; Lee & Weinberg 1977; Steigman & Turner 1985) and modified Newtonian dynamics (MOND; Milgrom 1983) were proposed as explanations of these observations. In MOND, the acceleration of test particles is modified, with respect to the Newtonian prediction, below a characteristic acceleration scale g₀ ∼ 10⁻¹⁰ m s⁻² to yield asymptotically constant speeds in rotation curves at large radii and low accelerations (Sanders & McGaugh 2002; Gentile et al. 2011; McGaugh et al. 2016; Lelli et al. 2017; Li et al. 2018), as observed. MOND also provides a correlation of this asymptotic speed with the total baryonic mass in the galaxy, that is, the baryonic Tully-Fisher relation (Tully & Fisher 1977; McGaugh et al. 2000). However, it has been argued that MOND cannot account for the entire missing mass observed in galaxy clusters (Sanders 2003), and today the more recent observations of merging clusters (Clowe et al. 2006) and the measurements of the cosmic microwave background (Skordis et al. 2006; Dodelson & Liguori 2006; Dodelson 2011) are considered by many as a challenge for MOND. A review of MOND and observations are available in Famaey & McGaugh (2012). Despite these known challenges, it is of obvious interest to investigate in detail the predictions of MOND for rotation curves beyond the asymptotic velocities at large radii.

Recently, the sample of galaxies in the SPARC database was compared to a MOND modified inertia (MI) model, and it was found that the fit residuals were Gaussian and of the expected size (McGaugh et al. 2016; Lelli et al. 2017; Li et al. 2018). However, in Frandsen & Petersen (2018) it was shown that data at small radii deviate significantly from the MOND MI fit. Since the sample of data at small radii is only a few hundred points, compared to the few thousand points at large radii, the discrepancy is only apparent in the residuals if these data are treated separately (see Fig. 5) or if galaxies are considered individually (Petersen & Frandsen 2020, see also Fig. 4). Investigating data points at small radii separately is well motivated as this is where the predictions of different models of the missing mass problem deviate significantly. In particular, MOND modified gravity (MG) and MOND MI models only deviate in their predictions in a definite manner at small radii, as we show here.

Early simulations of structure formation with cold collisionless DM particles and no baryons found universal cuspy Navarro-Frenk-White-like DM density profiles in halos ranging from dwarf galaxies to galaxy clusters (Navarro et al. 1996). This profile fits rotation curve data at large radii in galaxies and clusters, but not in all cases at small radii. The inferred DM densities from some observed clusters and gas-rich halo-dominated dwarf spirals is less steep (i.e., more cored) than the NFW profile in the inner regions (Flores & Primack 1994). This has become known as the cusp-core problem for DM. More recent DM-only simulations have found some systematic departures from the NFW profile and some diversity in the resulting rotation curves (Navarro et al. 2010). However, these DM-only simulations still show little variation in rotation curve profiles with the same asymptotic maximal rotational velocity (Oman et al. 2015), while observed rotation curves of dwarf galaxies do show such a variation. This has been termed the diversity problem. Whether the cusp-core problem or the diversity problem is a problem of DM, or rather of simulations with limited resolution and without the inclusion of baryons, remains a topic of debate since simulations that include baryonic feedback do find cored profiles (Read & Gilmore 2005; Teyssier et al. 2013; Di Cintio et al. 2014).

In this paper we identify a different related cusp-core challenge pertaining to MOND that is essentially the opposite of that for DM. To do so, we first provide a definition and classification of cuspy, neutral, and cored galaxies. The definition is based on acceleration curve geometry in g2 space following Frandsen & Petersen (2018) (i.e., the space of total centripetal accelerations, g_obs, versus the Newtonian centripetal acceleration from baryonic matter, g_bar) rather than on inferred DM density profiles. This classification is directly applicable to both MOND and DM models. Throughout this paper we refer to this new definition of cuspy and cored galaxies unless otherwise stated.

We show that MOND MG models, in the Bekenstein-Milgrom formulation, lead to cored acceleration curve geometries in g2 space and corresponding rotation speed curves – a consequence of the so-called solenoidal acceleration field in these models. We first illustrate this using analytical approximations (Brada & Milgrom 1994) and by investigating the curl of the solenoidal field, before we explicitly solve the MOND modified Poisson equation using the N-Mody code (Ciotti et al. 2006). In contrast, MOND MI curves provide a definition of ‘neutral’ curves, neither cuspy nor cored, as a benchmark.

A way to test MOND MG is therefore to look for cuspy galaxies according to our definition as well as cored galaxies that deviate from the specific cored geometries predicted by MOND MG models. We identify such galaxies in the SPARC database, extending our previous analyses of MOND MI models (Petersen & Frandsen 2020; Frandsen & Petersen 2018).

The paper is organized as follows. In Sect. 2 we present our definition and classification of cuspy, neutral, and cored galaxies based on acceleration curve geometry in g2 space. In Sect. 3 we show that MOND MG yields cored rotation curve geometries, which tend to neutral geometries for spherical mass distributions. We show this using general properties of the MOND solenoidal field, analytic approximations (Brada & Milgrom 1994), and full numerical solutions using the N-Mody code Ciotti et al. (2006). On the other hand, MOND MI universally leads to neutral rotation curve geometries, independent of the matter distribution. We note that Banik et al. (2018) also compares aspects of the analytical approximations in Brada & Milgrom (1994) to numerical simulations of the quasi-linear version of MOND MG, QUMOND. In Sect. 4 we identify cored and cuspy acceleration curve examples from the SPARC database of galaxies and compare them to the model predictions before giving our summary.

2. Classification of cuspy, cored, and neutral geometries in g2 space

We begin by reviewing the geometric characteristics of galactic centripetal acceleration curves, which we refer to as g2-space curves – the space of total predicted¹ centripetal accelerations, g_tot, versus the Newtonian acceleration from baryons, g_N – following Frandsen & Petersen (2018). Examples of MOND MG and MOND MI curves in g2 space are shown in the top panel of Fig. 1.

Fig. 1. Top panel: centripetal acceleration curves in g2 space with the quantities rtot, rN, and 𝒞±, used for classification in Table 1, shown. The solid grey line shows a MOND MI curve with the radius of maximum baryonic acceleration and maximum total acceleration rN = rtot indicated. Also, the curve segments 𝒞+ and 𝒞− coincide, so the curve area is 𝒜(𝒞) = 0. The dotted grey and dotted black curves show the 𝒞± curve segments of a MOND MG curve, using the Brada-Milgrom approximation in Eq. (9), with rtot >  rN and 𝒜(𝒞) > 0. For both, the baryonic matter is an infinitely thin exponential disc, Σ(r) = Σ0e−r/rd. The arrow indicates the direction of increasing radius along the curve. Bottom panel: corresponding rotation curves.

We first defined the radii r_tot and r_N as the locations at which the centripetal accelerations, g_tot and g_N, are maximum:

$$\begin{array}{c}\hfill g_{\mathrm{N}}(r_{\mathrm{N}})=\max \{g_{\mathrm{N}}(r)\},g_{\mathrm{tot}}(r_{\mathrm{tot}})=\max \{g_{\mathrm{tot}}(r)\}.\end{array}$$(1)

These radii are indicated in the top panel of Fig. 1 and coincide for MOND MI, while r_tot >  r_N for MOND MG. From these radii we defined the acceleration ratios,

$$\begin{array}{c}\hfill {\widehat{g}}_{\mathrm{tot}}(r)\equiv g_{\mathrm{tot}}(r)/g_{\mathrm{tot}}(r_{\mathrm{N}}),{\widehat{g}}_{\mathrm{N}}(r)\equiv g_{\mathrm{N}}(r)/g_{\mathrm{N}}(r_{\mathrm{N}}),\end{array}$$(2)

relevant for reducing systematic uncertainties when we study SPARC data (Frandsen & Petersen 2018).

We classified galaxies as cuspy if r_tot <  r_N, neutral if r_tot = r_N, and cored if r_tot >  r_N. More generally, we are interested in the relative location of the entire curve segment C⁺ at large radii and that at small radii C⁻ defined with respect to r_N:

$$\begin{array}{c}\hfill {\mathcal{C}}^{+}=\{(g_{\mathrm{N}}(r),g_{\mathrm{tot}}(r));r>r_{\mathrm{N}}\},{\mathcal{C}}^{-}=\{(g_{\mathrm{N}}(r),g_{\mathrm{tot}}(r));r<r_{\mathrm{N}}\}.\end{array}$$(3)

Galaxies are cored if 𝒞⁺ >  𝒞⁻, in the sense that the former curve segment lies above the latter, as seen for MOND MG in the top panel of Fig. 1; they are neutral if 𝒞⁺ = 𝒞⁻, as for the MOND MI in the top panel, and cuspy if 𝒞⁺ <  𝒞⁻. Using the quantities r_N, tot, C^± and the signed curve area 𝒜(𝒞), we summarize the definition of cuspy, neutral, and cored geometries in the first three rows of Table 1. The last three rows give the same using our notation for the same quantities in the data. Throughout this paper we refer to this new definition of cuspy, cored, and neutral galaxies unless otherwise stated.

Our definition is more general than that normally used for DM profiles, but a DM model with an NFW-like profile is also cuspy according to our definition, while that of DM with a quasi-isothermal profile is cored, as illustrated in Frandsen & Petersen (2018). From the right panel of Fig. 1, it is also seen that the MOND MG rotation speed curve indeed has a shallower approach to zero radius relative to the MOND MI rotation speed curve. This would correspond to a more cored density profile in the former case if it arose from DM. It is also seen that while the difference in g2-space curves is very significant, the effect is modest at the level of the rotation speed curve.

3. Cored and neutral MOND models in g2 space

In this section we show that MOND MG leads to cored geometries with r_tot >  r_N, 𝒞⁺ >  𝒞⁻ for isolated galaxies with axisymmetric mass distributions, while MOND MI, which universally leads to the geometries that we here term neutral, is characterized by r_tot = r_N, 𝒞⁺ = 𝒞⁻, as discussed in Frandsen & Petersen (2018). The resulting rotation speed curves of MOND MG are shallower than those of MOND MI, as exemplified in the bottom panel of Fig. 1. The cored MOND MG geometries reduce to the neutral MOND MI geometries for spherical mass distributions. To show the cored geometries of MOND MG, we start from analytical approximations and the simpler quasi-linear MOND MG models (Milgrom 2010) and finally solve the MOND MG curves explicitly using the N-MODY code Ciotti et al. (2006).

3.1. MOND modified gravity

In the Bekenstein-Milgrom formulation of MOND MG models (Bekenstein & Milgrom 1984), the total centripetal acceleration in a galaxy is determined via a modified Poisson equation for the MOND potential field, ψ,

$$\begin{array}{c}\hfill \mathbf{\nabla }·p(\mu (\frac{|\mathbf{\nabla }\psi |}{g_0})\mathbf{\nabla }\psi )=4\pi G\rho .\end{array}$$(4)

In this equation, −∇ψ = g_M is the MOND acceleration and g₀ ∼ 10⁻¹⁰ m s⁻² is a characteristic acceleration scale such that the interpolation function, μ(x), smoothly interpolates between the Newtonian regime, μ(x)≃1, for x ≫ 1 and the deep MONDian regime, μ(x)≃x, for x ≪ 1, with xμ(x) monotonic. The limiting behaviour of μ(x) for x ≫ 1 is clearly required to recover Newtonian dynamics, and the limiting behaviour in the deep MONDian limit leads to constant rotation curve speeds at large radii.

Using the Poisson equation for the Newtonian potential (Φ), 4πGρ = ∇²Φ = −∇ ⋅ g_N, where g_N is the acceleration predicted by Newtonian dynamics, the modified Poisson equation can be rewritten in terms of accelerations as

$$\begin{array}{c}\hfill \mu (\frac{g_M}{g_0}){\mathit{g}}_M={\mathit{g}}_{\mathrm{N}}+\mathit{S}=\mathit{q},{\mathit{g}}_{\mathrm{M}}=\nu (\frac{q}{g_0})\mathit{q},\end{array}$$(5)

where the solenoidal field, S, has zero divergence, ∇ ⋅ S = 0. The inverse interpolation function, ν(y), is defined such that ν(y)≡I⁻¹(y)/y with I(x) = xμ(x) = y. A number of interpolation functions, μ(x), and inverse interpolation functions, ν(y), have been considered in the literature (e.g., Begeman et al. 1991; Bekenstein 2004). Here, we consider two inverse interpolation functions: the ν₁(y) function used in the N-MODY code and the ν₂(y) proposed in Milgrom & Sanders (2008), McGaugh (2008), and Famaey & McGaugh (2012) and used to fit the SPARC galaxy data in McGaugh et al. (2016) and Lelli et al. (2017):

$$\begin{array}{c}\hfill \nu {(y)}_1=\sqrt{\frac{1}{2}(1+\sqrt{1+\frac{4}{y^2}})},\nu {(y)}_2=\frac{1}{1-e^{-\sqrt{y}}}.\end{array}$$(6)

The corresponding μ₁ function is ${\mu }_1(x)=\frac{x}{\sqrt{1+x^2}}$, and μ₂ is transcendental.

3.2. MOND modified inertia

For spherical matter distributions, S = 0 (Bekenstein & Milgrom 1984), and it follows that in this special case the Newtonian and MOND (g_MI) accelerations are related as

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{N}}=\mu (\frac{g_{\mathrm{MI}}}{g_0}){\mathit{g}}_{\mathrm{MI}};{\mathit{g}}_{\mathrm{MI}}=\nu (\frac{g_{\mathrm{N}}}{g_0}){\mathit{g}}_{\mathrm{N}}.\end{array}$$(7)

These relations with S = 0 also hold in so-called MOND MI models (Milgrom 1983) for any matter distributions. We therefore denote the MOND acceleration as g_MI in this case. The acceleration g_MI is generally not derivable from a potential as g_M defined via Eq. (4), but MOND MI is often used for comparing MOND with rotation curve data irrespective of the distinction between MOND MI and MOND MG. However, as seen above, the g2-space geometry of the two are distinctly different.

Since the function xμ(x) is monotonic, the function g_MI(g_N) is one-to-one. Then, since the Newtonian acceleration, g_N(r), goes to zero for both large and small radii, the g2-space curves, 𝒞, of MOND MI are closed curves with zero area. They are universally neutral geometries (i.e., independent of the underlying baryonic matter distribution and independent of the details of the interpolation function) according to the classification of Table 1 ².

3.3. The deep MONDian regimes at large and small radii

At large radii the solenoidal field, S, in Eq. (5) vanishes faster than the Newtonian acceleration, g_N, and so can be neglected. More precisely, the unit-less Newtonian acceleration y = g_N/g₀ falls off as y ∼ 1/r², while the solenoidal field, S, vanishes faster than 1/r³, as shown in Bekenstein & Milgrom (1984). Since the inverse interpolation function in the deep MONDian limit y → 0 behaves as ν(y)→y^−1/2, the large radius limit of the MOND acceleration is related to the Newtonian as

$$\begin{array}{c}\hfill g_{\mathrm{M}}\to \sqrt{g_{\mathrm{N}}{g}_0},\mathrm{for}r\to \infty .\end{array}$$(8)

The g2-space curves and rotation speed curves of MOND MI and MOND MG therefore always coincide at large radii. This is exemplified in the top panel of Fig. 1, where all curves coincide at large radii, corresponding to the g_N → 0 parts of the 𝒞⁺ curve segments (the upper curve segments).

3.3.1. Analytical approximations

To study the deep MONDian limit y → 0 at small radii, it is instructive to consider an analytic approximation to MOND MG for infinitely thin discs. In particular, we considered an exponential disc with a surface mass density Σ(r) = Σ₀e^−r/r_d. In this case, an approximate expression for the resulting centripetal acceleration g_BM, r in MOND MG was given by Brada and Milgrom in Brada & Milgrom (1994). Taking Eq. (7) as an approximation for the MOND MG acceleration outside the disc and then taking into account the discontinuity of the z component of the acceleration inside the disc, one finds the radial acceleration in the plane of the infinitely thin disc:

$$\begin{array}{c}\hfill g_{\mathrm{BM},r}(g_{\mathrm{N}},r)=\nu \left(\frac{g_{\mathrm{N}}^{+}}{g_0}\right)g_{\mathrm{N},r},g_{{\mathrm{N}}^{+}}=\sqrt{g_{\mathrm{N},r}^2+{(2\pi G\mathrm{\Sigma }(r))}^2},\end{array}$$(9)

where g_N⁺ is the total acceleration in the disc, which includes the discontinuity in the z component by taking the limit z → 0⁺ from the upper half of the plane. Given the limiting behaviour of the density Σ(r) at large and small radii, we find the centripetal acceleration in the two deep MONDian regimes to be

$$\begin{array}{cc}\hfill g_{\mathrm{BM},\mathrm{r}}& =\nu \left(\frac{g_{\mathrm{N}}^{+}}{g_0}\right)g_{\mathrm{N},r}\to \sqrt{g_{\mathrm{N},r}{g}_0},\mathrm{for}r\to \infty ,\hfill \\ \hfill g_{\mathrm{BM},\mathrm{r}}& =\nu \left(\frac{g_{\mathrm{N}}^{+}}{g_0}\right)g_{\mathrm{N},r}\to \sqrt{\frac{g_0}{{\mathrm{\Sigma }}_0}}{g}_{\mathrm{N},r},\mathrm{for}r\to 0.\hfill \end{array}$$(10)

At large radii, this MOND MG approximation coincides with MOND MI, as it should, while at small radii it is reduced by a factor $\sqrt{\frac{g_{\mathrm{N},r}}{{\mathrm{\Sigma }}_0}}<1$ relative to MOND MI. In the deep MONDian regime at small radii, the MOND and Newtonian accelerations are linearly related, as opposed to the square root relation in the deep MONDian regime at large radii. We show this MOND MG approximation in Fig. 1 as the dotted acceleration (left panel) and speed curves (right panel) and in Fig. 3 as the dotted black curve. The square root and linear behaviour of g_BM, r as a function of g_N, r in the deep MONDian regimes at small and large radii, respectively, are clearly seen in the figures. It follows that, in this approximation, MOND MG curves are cored with r_tot >  r_N, 𝒞⁺ >  𝒞⁻, and 𝒜(𝒞) > 0.

It is also possible to start from a spherical potential for which the MOND MG acceleration coincides with the MOND MI exactly. By adding an axisymmetric perturbation, one can compute the solenoidal field, S, as a function of this perturbation. This is done in Ciotti et al. (2006). We now discuss how the cored geometry of MOND MG arises beyond the infinitely thin disc approximation or the approximation of a nearly spherical mass distribution.

3.3.2. Quasi-linear MOND modified gravity

Before solving the full MOND MG geometries, it is also instructive to consider QUMOND, the quasi-linear version of MOND MG (Milgrom 2010), to see the origin of the cored geometry. The QUMOND acceleration, g_QM, is obtained by starting from the (pristine) MOND MI acceleration, g_MI, in the right-hand side of Eq. (7) and then adding a solenoidal field, σ = ∇ × A, to enforce that the resulting QUMOND acceleration, g_QM, has zero curl. This ensures that the QUMOND acceleration is derivable from a potential:

$$\begin{array}{c}\hfill {\mathit{g}}_{\mathrm{QM}}=\nu \left(\frac{g_{\mathrm{N}}}{g_0}\right){\mathit{g}}_{\mathrm{N}}+\mathit{\sigma },\mathrm{\nabla }\times \mathit{\sigma }=-\frac{{\nu }^{\prime }}{g_0}\mathrm{\nabla }{g}_{\mathrm{N}}\times {\mathit{g}}_{\mathrm{N}},\end{array}$$(11)

where the curl of σ follows from the requirement ∇ × g_QM = 0 such that g_QM = −∇ψ_QM. This allows us to determine ∇ × σ straightforwardly in terms of the Newtonian potential.

Due to the axisymmetry of the matter distribution, the curl ∇ × σ is purely in the azimuthal direction, and, using the second identity in Eq. (11), we can also determine the sign of ∇ × σ in the azimuthal direction to be negative. This is shown in Fig. 2 for a Miyamoto-Nagai (MN) potential, which is given in Eq. (13) and reviewed in Sect. 3.4. The scale length (a), scale height (b), and mass scale (M) of the MN model used in the figure are a = 1, b = 0.5, and M = 1 (bottom left). The negative sign arises as follows: For z ≥ 0, the gradient ∇g_N is dominantly in the z direction for a non-spherical, axisymmetric mass distribution (it is zero in the plane by symmetry unless the disc is infinitely thin), while g_N is dominantly in the negative r direction. This is illustrated for the MN model in Fig. 2 in the left and middle panels, respectively. It is particularly clear at the line of maximum radial acceleration, ∂_rg_N = 0, near r = 1.1, where ∇g_N points entirely in the positive z direction and g_N is dominantly in the negative radial direction. Then, since ν′ < 0 in the MONDian regimes at small and large radii, we have that ∇ × σ is in the negative azimuthal direction, as seen in the right-hand panel of the figure, and, consequently, σ circulates in the anticlockwise direction in the (r, z) plane.

Fig. 2. Gradient of the norm of the Newtonian acceleration, ∇gN (left panel), the Newtonian acceleration vector, gN (middle panel), and the curl of the solenoidal field, ∇ × σ (right panel), in QUMOND for an MN disc model with mass parameter MMN = 1, scale height b = 0.5, and scale length a = 1. From the sign of ∇ × σ, we infer that the centripetal component of the curl acceleration is in the opposite direction of the Newtonian acceleration vector, gN, in the disc near z = 0.

At large radii σ is negligible, but at small radii the radial component of g_QM is reduced as compared to the radial component of g_MI because the radial acceleration from the solenoidal field, σ, is in the opposite direction to the Newtonian acceleration, g_N. It follows, using the classification of Table 1, that the QUMOND MG geometries are cored.

3.3.3. Solenoidal field in MOND modified gravity

The arguments above for QUMOND may also be applied to the solenoidal field, S, of the Bekenstein-Milgrom MOND MG (note the opposite sign convention!). The curl ∇ × S may be expressed from either of the two identities in Eq. (12), using ∇ × g_N = ∇ × g_M = 0, as

$$\begin{array}{c}\hfill \mathrm{\nabla }\times \mathit{S}={\mu }^{\prime }\mathrm{\nabla }{g}_M\times {\mathit{g}}_M=-\frac{{\nu }^{\prime }}{\nu g_0}\mathrm{\nabla }q\times \mathit{q}.\end{array}$$(12)

Due to the axisymmetry, the curl, ∇ × S, is again purely in the azimuthal direction. We can argue that the sign of ∇ × S is in the negative azimuthal direction such that S circulates in the anticlockwise direction by repeating the arguments for QUMOND: At large radii, the solenoidal field is negligible, and as we move towards smaller radii we approximate q ≃ g_N in the above formula and repeat the arguments above for QUMOND (i.e. near the z = 0 plane the radial acceleration is reduced as compared to that in MOND MI). It follows that the MOND MG geometries are cored, as in QUMOND, using the classification of Table 1. Example numerical solutions of Eq. (4), using the N-MODY code, are shown in the top-left panel of Fig. 3 (coloured dotted curves) as a function of the parameter b, which controls the departure from the spherical limit. We now discuss these numerical solutions.

Fig. 3. Top-left panel: centripetal acceleration curves in g2 space of MOND MI (solid grey curve) and MOND MG in the Brada-Milgrom approximation (black dotted curve) with acceleration scale g0 = 1.2 × 10−10 m s−2 for a 3MN model of an infinitely thin exponential disc galaxy with M = 1.2 × 1010 M⊙ and a scale length of rd = 3.5 kpc. The 3MN parameters are given in the appendix. The dotted coloured curves are the full MOND MG solution of the same model but with the scale height parameter of the 3MN model varied: b = 3 (green), b = 1 (red), and b = 0.3 (blue) (b → 0 is the disc limit, and b → ∞ is the spherical limit). Finally, the solid orange line gives the purely Newtonian acceleration. Top-right panel: corresponding rotation curves showing the full MOND MG solution only for b = 0.3. Bottom-left panel: corresponding curl field S for the b = 3 case, showing how the curl field at z = 0 is in the opposite direction to the Newtonian one and how it reduces the radial MOND MG acceleration compared to the MOND MI. Bottom-right panel: difference in curl fields for the most disc-like (b = 0.3) and most spherical (b = 3) parameter values, showing that the solenoidal field reduces the radial MOND MG acceleration in the z = 0 plane most for the most disc-like mass distributions.

3.4. Numerical solutions of MOND modified gravity

In order to study numerical solutions of the full MOND MG geometries in g2 space, it is useful to consider MN potentials, which interpolate between axisymmetric potentials and spherical potentials. The MN potential, Φ_MN, is given by

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{MN}}(r,z)=\frac{-G{M}_{\mathrm{MN}}}{\sqrt{r^2+{(a+\sqrt{z^2+b^2})}^2}},\end{array}$$(13)

where M_MN is a mass parameter, b is the vertical scale height, and a is the radial scale length (Miyamoto & Nagai 1975). In the limit a = 0, the potential is spherically symmetric and reduces to the Plummer potential. For a ≠ 0 the potential is axisymmetric (and symmetric about the z = 0 plane), and in the limit b = 0 the potential reduces to the potential of an infinitely thin Kuzmin disc, for which the MOND MG solution is known exactly.

However, instead of a single MN model, we used a sum of three Miyamoto-Nagai potentials (3MN models), which can be used to approximate the potential of both thick and thin exponential discs (Smith et al. 2015), with densities of the form

$$\begin{array}{c}\hfill \rho (r,z)={\rho }_{0}exp(-r/r_d)exp(-|z|/z_0),{\rho }_0=\frac{M}{4\pi h_z{r}_d^2},\end{array}$$(14)

where r_d is the radial scale length and z₀ is the vertical scale length. This allows us to compare numerical solutions of MOND MG to both the Brada-Milgrom approximation for infinitely thin exponential discs and MOND MI in the spherical limit. The 3MN potentials we used are therefore of the form:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{3\mathrm{MN}}(r,z)={\displaystyle \sum _{i=1}^3}\frac{-G{M}_{\mathrm{MN},i}}{\sqrt{r^2+{(a_i+\sqrt{z^2+b_i^2})}^2}}.\end{array}$$(15)

We discuss details of the 3MN approximation in the appendix.

As an example, we started from an infinitely thin exponential disc galaxy with a mass of M = 1.2 × 10¹⁰ M_⊙ and a scale length of r_d = 3.5 kpc and took g₀ = 1.2 × 10⁻¹⁰ m s⁻². The corresponding 3MN model parameters, M_3MN, b = 0, a_i, M_i, are given in the appendix and are used to compute the MOND MI g_MI curve (solid grey line) and MOND MG in the Brada-Milgrom approximation, g_BM (dotted black line), in the top-left panel of Fig. 3. By dialling the parameter b up, we can make the mass distribution more spherical, with a/b → 0 the spherical limit (corresponding to increasing the scale height, z₀, in the thick exponential disc above).

The coloured dotted curves show the full solution of the MOND MG curve, using the static solver of the N-MODY code for the 3MN models with increasing values of b from b = 0.3 (blue), to b = 1 (orange), to b = 3 (green).

The figure demonstrates the properties discussed previously: The MOND MI and MOND MG curves coincide at large radii, but the MOND MG curves are more cored at small radii and the core grows in proportion to the departure from spherical symmetry (as measured by the ratio a/b here) towards the Brada-Milgrom approximation in the thin disc limit a/b → ∞. We show the rotation speed curves corresponding the MOND MI and Brada-Milgrom curves along with the full numerical solution for b = 0.3 in the top-right panel. Finally, it is seen how the approach to zero acceleration of the curves follows the radial limiting behaviour given in Eq. (10). We also note that increasing the mass of the disc, such that the accelerations reach the Newtonian regime, simply means the curves are stretched and the core effect at small radii is less pronounced.

In the bottom-left panel we show the corresponding solenoidal field, S, for the b = 3 example. The curl field indeed circulates in the anticlockwise direction. Finally, in the bottom-right panel of the figure, we show that the difference between the b = 0.3 and b = 3 curl fields is positive near z = 0; therefore, the least spherical potential leads to the biggest curl field along the z = 0 plane and thus the biggest core effect. We can therefore take the Brada-Milgrom approximation as a good approximation of MOND MG for thin discs. For more spherical mass distributions, the Brada-Milgrom will overestimate the core, so for a given galaxy we can take the Brada-Milgrom formula and apply it to a thin disc of equivalent mass to get an upper limit on the core produced by MOND MG in the galaxy.

In summary, we have shown that MOND MG leads to cored g2-space geometries (and corresponding rotation curves) for isolated galaxies and that the core effect is then smallest for the most spherical systems, for example dwarf spheroidals. MOND MI, which is also the limit of MOND MG for fully spherical baryonic mass distributions, leads to neutral geometries.

4. MOND models and SPARC galaxy data in g2 space

We now compare MOND model acceleration curves to data from the SPARC galaxy database (Lelli et al. 2016). We focus on galaxies that exhibit cuspy and cored geometries, respectively.

The geometry of MOND MG is cored relative to MOND MI. Therefore, the (simpler) fits of MOND MI to cuspy galaxies provide an upper limit on the goodness of fit for MOND MG to such galaxies. Conversely from Fig. 3, we can use the Brada-Milgrom approximation as an upper limit on how cored MOND MG can be (i.e., how big the enclosed area 𝒜(C) is). Our data analysis below follows that in Frandsen & Petersen (2018) and is detailed in Appendix A.

4.1. Cored and cuspy SPARC galaxies

In the top-left panel of Fig. 4, we show galaxies from SPARC with cored g2-space geometries (colour legend in the figure identifies the individual galaxies) along with the MI model with the interpolation functions in Eq. (6) and with the best-fit value g₀ = 1.2 × 10⁻¹⁰ m s⁻² from Lelli et al. (2016) (solid and dashed grey lines). The galaxies were selected as cored only by requiring r_obs, G >  r_bar, G and by having a large ${\chi }_G^2$ value for the fit of data from each galaxy to the MOND MI curve in Eq. (6). Below, we discuss the fits in more detail. In the top-right panel we show the same for galaxies with cuspy g2-space geometries, selected as cuspy only by requiring r_obs, G < r_bar, G and, again, by having a large ${\chi }_G^2$ for the fit of the data from each galaxy to the curve in Eq. (6).

Fig. 4. Top-left panel: acceleration curves of galaxies from SPARC in g2 space with cored geometry robs >  rbar. Top-right panel: same but for cuspy galaxies with robs <  rbar. Also shown in both panels are the MOND MI curves for the two interpolation functions considered in Eq. (6). Middle-left panel: corresponding model curves for the cored galaxies from MOND MG curves using the Brada-Milgrom approximation with the gbar and surface density Σ(r) values from SPARC. Middle-right panel: same as left but for the cuspy galaxies. For reference, the MOND MI curve is also shown again. Bottom-left panel: corresponding rotation speed curve data for two of the cored galaxies, compared to the MOND MI and MOND MG curves. Bottom-right panel: same as left but for two cuspy galaxies

In the middle row of panels we show the predicted MOND MG curves of the same cored and cuspy galaxies along with the MOND MI curves for reference. We obtained the MOND MG curves using the Brada-Milgrom approximation in Eq. (9), with g_N replaced by the data g_bar(r_j, G) and with Σ(r) replaced by the disc surface densities Σ(r_j, G) from SPARC.

It is visually clear that the MOND MG curves yield poor fits to both the cored and cuspy galaxies, as we quantify further below. For the cuspy galaxies the MOND MG is, by nature of its cored geometry, a worse description than MOND MI. For the cored galaxies, the MOND MG curves are everywhere too close to the MOND MI limit to describe the cored geometry of the data. In particular, the MOND MG curves are more cored for the cuspy galaxies than they are for the cored ones because the latter have larger stellar surface mass densities.

In the bottom row of panels we show the corresponding rotation speed curves for two of the cored and two of the cuspy galaxies compared to the predictions from MOND MI and the Brada-Milgrom approximation of MOND MG.

In Fig. 5 we show the distribution of normalized residuals, R(r_j, G), from Eq. (A.3) of the SPARC data after quality cuts and the δv_obs/v_obs <  0.1 cut employed in McGaugh et al. (2016), Lelli et al. (2017), and Lelli et al. (2016) with respect to the MOND MI prediction with g₀ = 1.2 × 10⁻¹⁰ m s⁻². The grey histograms, highlighted with black dots in the middle of the bins, show the residuals of all SPARC data points at large radii (i.e., with r_j, G >  r_bar, G). These points are seen to follow the Gaussian of unit width superimposed on the figure, as expected of a good fit. The left panel also shows the residuals of data at small radii, meaning r_j, G ≤ r_bar, G from the (cored) galaxies with r_obs, G >  r_bar, G (red histogram). These residuals are not Gaussian but skewed towards large negative residuals, which is consistent with these points being below the MOND MI prediction in general, as is the case for our example galaxies. This is also what MOND MG would in general predict they should be; however, for our examples in Fig. 4, the specific cores of MOND MG do not match the data well. In the right panel we also show the residuals of data at small radii, r_j, G ≤ r_bar, G, from (cuspy) galaxies with r_obs, G <  r_bar, G (blue histogram). These residuals are skewed towards large positive residuals, consistent with these points lying above the MOND MI predictions and therefore also above the MOND MG predictions in g2 space.

Fig. 5. Top panel: distribution of the normalized residuals R(rj, G) in Eq. (A.3) of SPARC data with respect to MOND MI with g0 = 1.2 × 10−10 m s−2. The grey histograms with black dots at each bin top shown on both panels are the 1983 SPARC data points at large radii, i.e., rj, G >  rbar, G, after imposing the basic quality criteria and the δvobs/vobs <  0.1 cut employed in McGaugh et al. (2016), Lelli et al. (2017), and Lelli et al. (2016). Also shown in the top panel are normalized residuals of data points at small radii, rj, G ≤ rbar, G, from the (cored) galaxies only with robs, G >  rbar, G (red histogram). Bottom panel: points at small radii, rj, G <  rbar, G, from (cuspy) galaxies only with robs, G <  rbar, G (cuspy; blue histogram).

We note that seven of the galaxies in SPARC with steeply rising (cuspy) rotation curves are starburst dwarf galaxies, for which data may not faithfully represent the underlying gravitational potential, as discussed in Read et al. (2016) and Santos-Santos et al. (2018). However, some of those galaxies are eliminated by the data quality cut employed here, and they are not among the galaxies presented in Fig. 4.

4.2. Model fits to data

To reduce systematic uncertainties in the data before performing quantitative fits of MOND to the galaxies in Fig. 4, we defined the radius of maximum baryonic acceleration in the data (in analogy with the model r_N above) and acceleration ratios as:

$$\begin{array}{cc}& g_{\mathrm{bar}}(r_{\mathrm{bar},G})=\max \{g_{\mathrm{bar}}(r_{j,G})\};{\widehat{g}}_{\mathrm{obs}}(r_{j,G})\equiv \frac{g_{\mathrm{obs}}(r_{j,G})}{g_{\mathrm{obs}}(r_{\mathrm{bar},G})},\hfill \\ \hfill & {\widehat{g}}_{\mathrm{bar}}(r_{j,G})\equiv \frac{g_{\mathrm{bar}}(r_{j,G})}{g_{\mathrm{bar}}(r_{\mathrm{bar},G})}.\hfill \end{array}$$(16)

These ratios, introduced in Frandsen & Petersen (2018), eliminate the systematic uncertainties on galactic distance and inclination angle for g_obs, and they significantly reduce the systematic error from gas measurements and mass to-light ratios for g_bar. The resulting uncertainties, $\delta {\widehat{g}}_{\mathrm{bar}}(r_{j,G})$ and $\delta {\widehat{g}}_{\mathrm{obs}}(r_{j,G})$, are given in Frandsen & Petersen (2018) and reproduced in Appendix B. Using the acceleration ratios, we constructed the ${\chi }_G^2$ of each galaxy compared to the model acceleration,

$$\begin{array}{c}\hfill {\chi }_G^2={\displaystyle \sum _{j,j^{\prime }}}({\widehat{g}}_{\mathrm{obs}}(r_{j,G})-{\widehat{g}}_M(r_{j,G}))V_{j{j}^{\prime }}^{-1}({\widehat{g}}_{\mathrm{obs}}(r_{j^{\prime },G})-{\widehat{g}}_M(r_{j^{\prime },G})),\end{array}$$(17)

where ${\widehat{g}}_M(r_{j,G})={\widehat{g}}_M({\widehat{g}}_{\mathrm{bar}}(r_{j,G}))$ is the MOND model prediction with the data point ${\widehat{g}}_{\mathrm{bar}}(r_{j,G})$ as input, and the inverse variance matrix, $V_{j{j}^{\prime }}^{-1}$, given in in Appendix B, takes into account the systematic uncertainty from the normalization point common to all acceleration ratios within a single galaxy. We neglected the uncertainties in ${\widehat{g}}_{\mathrm{bar}}(r_{\mathrm{bar},G}),$ which are small compared to those in ${\widehat{g}}_{\mathrm{obs}}(r_{j,G})$ for most data points, as seen in Figs. 7 and 6, where the acceleration ratios with errors are shown.

Fig. 6. Acceleration ratio curves $({\widehat{g}}_{\mathrm{bar}},{\widehat{g}}_{\mathrm{obs}})$ (points with errors) in $\widehat{g}2$ space for three of the cuspy SPARC galaxies shown in the right-hand panels of Fig. 4, compared to the MOND MI curve with g0 = 1.2 × 10−10 (solid grey), the MOND MI curve with the best-fit value of g0, min for each galaxy (dashed grey), and MOND MG in the Brada-Milgrom approximation curves (points without errors). Also shown are the χ2/d.o.f. values of the MOND MI curve with g0 = 1.2 × 10−10.

Fig. 7. $\widehat{g}2$-space data points (points with errors) of three of the cored SPARC galaxies shown in the left-hand panel of Fig. 4, compared to the MOND MI curve with g0 = 1.2 × 10−10 (solid grey), the MOND MI curve with the best-fit value of g0 for each galaxy (dashed grey), and MOND MG (Brada-Milgrom approximation) curves (points without errors). We note that most ${\widehat{g}}_{\mathrm{bar}}$ uncertainties are small for most cases because the ratios are constructed to reduce systematic uncertainties. Also shown are the χ2/d.o.f. values for the dashed best-fit MOND MI curves.

4.2.1. Cuspy SPARC galaxies

We first studied three of the cuspy galaxies shown in the top-right panel of Fig. 4. Their corresponding data curves – with errors in the normalized $\widehat{g}$ variables and after imposing the data cut δv_obs/v_obs <  0.1 used in Lelli et al. (2017) and Lelli et al. (2016) – are shown in Fig. 6. Also shown are the MOND MI curves with g₀ = 1.2 × 10⁻¹⁰ m s⁻² (solid grey) and with the best-fit value g_0, min that minimizes the ${\chi }_G^2$ in Eq. (17) (dashed grey). Finally, the MOND MG model curve using the Brada-Milgrom approximation is shown as points without errors in the same colour as the data.

The uncertainties on ${\widehat{g}}_{\mathrm{bar}}$ are indeed small compared to those on ${\widehat{g}}_{\mathrm{obs}}$ for most points. We could therefore compute the ${\chi }_G^2$ value for each of these galaxies with respect to MOND using Eq. (17). We give the ${\chi }_G^2$ value of MOND MI with g₀ = 1.2 × 10⁻¹⁰ m s⁻² fixed (McGaugh et al. 2016) in the second row of Table 2 and the minimum ${\chi }_G^2$ with the corresponding g_0, min value in the third row. As already discussed, and as is clear from the figures, the ${\chi }_G^2$ values of MOND MG will be larger or equal to those of MOND MI quoted here.

4.2.2. Cored SPARC galaxies

We next considered fits of MOND to three of the cored galaxies shown in g2 space in the top-left panel of Fig. 4. The corresponding acceleration ratio curves are shown in Fig. 7. The ${\chi }_G^2$ values of these galaxies with respect to MOND MI, analogous to those above for cuspy galaxies, are listed in Table 3.

Again, we find large ${\chi }_G^2$ values for the fits to MOND MI, and, given the very little difference between this and MOND MG seen in the figures, we can also take these values for MOND MG.

As summarized in Tables 2 and 3, the MOND fits to individual cuspy and cored galaxies return large χ² values and a significant variation in g₀. These best-fit values of g₀ deviate considerably from the value g₀ = 1.2 × 10⁻¹⁰ m s⁻², which was found in Lelli et al. (2017) as a best fit to the entire SPARC dataset.

These deviations can be understood following the discussion in Petersen & Frandsen (2020): At large radii, MOND MI and MOND MG are indistinguishable, as shown in Fig. 3. (This is so since MOND MG reduces to MI in the spherical limit, which is approximately valid at large distances). They both predict $g_{\mathrm{tot}}\simeq \sqrt{g_{\mathrm{N}}}$ at these radii, and we should expect the models to reproduce data well by virtue of their original construction. By contrast, MOND MI and MOND MG differ distinctly at radii r ≲ r_bar, where the difference between cuspy and cored acceleration curve geometries particularly manifests itself (it is also where the difference between cuspy and cored DM density profiles manifests itself).

The SPARC database is dominated by points at large radii, r >  r_bar (the 1983 points in grey in Fig. 5), which are indeed in good agreement with both MOND MI and MOND MG. Compared to MG, there are 385 points at small radii, r ≤ r_bar (the red and blue points shown in Fig. 5 ), which are in poor agreement with MOND; Fig. 5 displays the large non-Gaussian residuals of these data points with respect to MOND MI with g₀ = 1.2 × 10⁻¹⁰ m s⁻².

Since most of the data points are at large radii, r >  r_bar, in SPARC, a fit to all these data points that does not retain radial information, as the one in Lelli et al. (2017), will also mostly test the behaviour of acceleration curves at large radii. However, radial information from each galaxy is needed in order to fully test any model and to distinguish between different models – for example, by simply analysing individual galaxies, as we do in this study.

5. Summary

In this paper we have provided and employed a new definition of cuspy and cored galactic acceleration curve geometries that is applicable to both DM and MOND models of missing mass. Cuspy and cored acceleration curves are defined relative to the curves of MOND MI, which we take as a neutral reference, and the classification in the space of baryonic and total accelerations (g2 space) is that proposed in Frandsen & Petersen (2018). This new definition is summarized in Table 1, illustrated in Fig. 1, and used throughout, unless otherwise stated.

Examples of cored and cuspy acceleration curves in g2 space, and corresponding rotation speed curves, from the SPARC database are shown in the top and bottom panels of Fig. 4. We fit MOND MI and MOND MG to some of these galaxies in Figs. 6 and 7, with the intrinsic MOND acceleration scale g₀, defined in Eq. (4), as fit parameter. To eliminate systematic uncertainties from galaxy distance and inclination angle and to reduce the systematic uncertainty from mass-to-light ratios, the fit is performed on ratios of accelerations ${\widehat{g}}_{\mathrm{bar}}$ and ${\widehat{g}}_{\mathrm{obs}}$ defined in Eq. (16). As summarized in Tables 2 and 3, the fits return large χ² values and significant variation in the best-fit values of g₀ between different galaxies. Also, these best-fit values of g₀ deviate considerably from the value g₀ = 1.2 × 10⁻¹⁰ m s⁻² that was found in Lelli et al. (2017) as a best fit to the entire SPARC dataset. These deviations are visually clear from Fig. 4 and are in line with our previous finding: that the neutral geometry of MOND MI, specifically the prediction g_tot(r_tot) = g_tot(r_N), is in tension with the full SPARC dataset of circa 150 galaxies Frandsen & Petersen (2018), independent of the MOND interpolation function chosen.

The cored acceleration curves from MOND MG models are most pronounced for galaxies with the most disc-like baryonic matter distribution, and, in the limit of spherical distributions, they reduce to the neutral curves of MOND MI with the functional relation g_tot = g_tot(g_N). This is shown in Fig. 3 and discussed in Sect. 3. The cored and neutral curves of MOND can be contrasted with, for example, the cuspy curves from NFW DM density profiles that arise in N-body simulations of DM structure formation without baryons.

Based on our results and our new definition of cuspy, cored, and neutral acceleration curve geometries, we have elucidated an observational challenge for MOND, which we refer to as the cusp-core challenge: MOND MG leads to cored acceleration curves for isolated galaxies as a consequence of the solenoidal field, S, arising from the modified Poisson equation given in Eq. (5). MOND MI leads to neutral acceleration curves for isolated galaxies as a consequence of the one-to-one functional relationship between the Newtonian and MOND MI accelerations given in Eq. (7). In particular, neither allows for cuspy acceleration curves. By contrast, as illustrated in the top panels of Fig. 4, galactic acceleration curves inferred from observations display a diversity of cored and neutral as well as cuspy acceleration curve geometries.

This cusp-core challenge we have identified for MOND is analogous to – but essentially the inverse of – the DM cusp-core problem (Flores & Primack 1994), in the sense that DM-only N-body simulations lead to NFW-like cuspy DM density profiles (Navarro et al. 1996). By contrast, a diversity of cored and cuspy density profiles are inferred from observations.

Baryonic feedback from supernovae can change cuspy DM profiles into cored ones in some cases (Read & Gilmore 2005; Teyssier et al. 2013; Di Cintio et al. 2014). In the future, it would be interesting to investigate whether, for example, the external field effect in MOND or a radial dependence of the mass-to-light ratios impact and potentially alleviate the problem posed by the cusp-core challenge for MOND.

Acknowledgments

We thank J. Read, W.-C. Huang and J. Smirnov for discussions and comments on the draft. The authors acknowledge partial funding from The Council For Independent Research, grant number DFF 6108-00623. The CP3-Origins center is partially funded by the Danish National Research Foundation, grant number DNRF90.

References

Appendix A: Galaxy data analysis

The SPARC database provides the observed rotational velocities, v_obs(r_j, G), as well as the inferred contributions to the baryonic velocity from the disc, bulge, and gas in the galaxy, v_disc(r_j, G), v_bul(r_j, G), and v_gas(r_j, G), with r_j, G the jth measured radius in the galaxy, G. The accelerations g_obs(r_j, G) and g_bar(r_j, G) are then

$$\begin{array}{cc}& g_{\mathrm{obs}}(r_{j,G})={\displaystyle \frac{v_{\mathrm{obs}}^2(r_{j,G})}{r_{j,G}}},\hfill \\ \hfill & g_{\mathrm{bar}}(r_{j,G})=\frac{|v_{\mathrm{gas}}(r_{j,G})|v_{\mathrm{gas}}(r_{j,G})}{r_{j,G}}+{\displaystyle \sum _{i=\mathrm{disc},\mathrm{bulge}}}\frac{{\widehat{\mathrm{{\rm Y}}}}_{\mathrm{i},\mathrm{G}}{v}_{\mathrm{i}}{(r_{j,G})}^2}{r_{j,G}},\hfill \end{array}$$(A.1)

where ${\widehat{\mathrm{{\rm Y}}}}_{\mathrm{disc},\mathrm{G}}$ and ${\widehat{\mathrm{{\rm Y}}}}_{\mathrm{bulge},\mathrm{G}}$ are unit-less mass-to-light ratio rescaling factors. They are included following McGaugh et al. (2016), Lelli et al. (2017), Lelli et al. (2016), and Li et al. (2018) because the speed data, v_i(r_j, G), listed in the SPARC database are all quoted assuming reference mass-to-light ratios equal to that of the Sun, Υ_i, G = M_⊙/L_⊙. Following McGaugh et al. (2016), Lelli et al. (2017), Lelli et al. (2016), and Li et al. (2018), we take the central values, ${\widehat{\mathrm{{\rm Y}}}}_{\mathrm{disc},\mathrm{G}}=0.5$ and ${\widehat{\mathrm{{\rm Y}}}}_{\mathrm{bulge},\mathrm{G}}=0.7$, including a 25% relative uncertainty on both, $\delta \widehat{\mathrm{{\rm Y}}}=0.25\widehat{\mathrm{{\rm Y}}}$.

We note the absolute value on v_gas. This appears because for some galaxies in the SPARC database v_gas is reconstructed as negative in the innermost regions of the galaxies. This happens when the gas distribution has a large central depression and the gravitational force from the material in the outer regions is stronger than the force from the material in the inner regions. In this case, the gas acts to reduce rather than enhance the total centripetal acceleration, and this is taken into account in the accelerations computed from SPARC data via negative v_gas values for these radii as discussed in Lelli et al. (2016). The galaxy UGC03580, shown in Figs. 4 and 6, and the galaxy UGC08286, shown in Fig. 4 only, has such negative v_gas values. All the other galaxies used here do not.

The SPARC database also provides the corresponding (random) uncertainties δv_obs(r_j, G), as well as the uncertainties δi_G and δD_G on the galaxy inclination angle, i_G, and distance, D_G. Following Lelli et al. (2017), we further adopted a 10% uncertainty on the HI flux calibration, which translated into a 10% uncertainty on the gas accelerations (i.e. δg_gas(r_j, G) = 0.1g_gas(r_j, G)). From these uncertainty contributions, we can compute the full δg_bar, δg_obs uncertainties:

$$\begin{array}{cc}\hfill \delta g_{\mathrm{obs}}(r_{j,G})& =g_{\mathrm{obs}}(r_{j,G})\sqrt{[{\displaystyle \frac{2\delta v_{\mathrm{obs}}(r_{j,G})}{v_{\mathrm{obs}}(r_{j,G})}}{]}^2+[{\displaystyle \frac{2\delta i_G}{tan(i_G)}}{]}^2+[{\displaystyle \frac{\delta D_G}{D_G}}{]}^2},\hfill \\ \hfill \delta g_{\mathrm{bar}}(r_{j,G})& =g_{\mathit{bar}}(r_{j,G})\sqrt{(\frac{\delta g_{\mathrm{gas}}(r_{j,G})}{g_{\mathrm{bar}}(r_{j,G})}{)}^2+{\displaystyle \sum _{k=disc,bulge}}(\frac{v_k^2(r_{j,G})\delta {\mathrm{{\rm Y}}}_{k,G}}{v_{\mathit{bar}}^2(r_{j,G})}{)}^2},\hfill \end{array}$$(A.2)

where we note that the inferred g_bar(r_j, G) are independent of distance, D_G, and inclination angle, i_G (Li et al. 2018). We treat the uncertainties δv_obs as random Gaussian-distributed uncertainties for each data point, while the remaining uncertainties, δi_G, δD_G, δΥ_{disc, bulge, G}, and δg_gas(r_j, G), are systematic within each galaxy, meaning they rescale all data points within a given galaxy in one direction.

The normalized residuals that enter into the χ² function, taking the errors in both the observed and baryonic accelerations into account, can be approximated as (Orear 1982)

$$\begin{array}{c}\hfill R(r_{j,G})=\frac{(g_{\mathrm{obs}}(r_{j,G})-g_{\mathit{MI}}(r_{j,G}))}{\sqrt{\delta g_{\mathrm{obs}}{(r_{j,G})}^2+g_{\mathit{MI}}^{\prime }{(r_{j,G})}^2{g}_{\mathrm{abs}}{(r_{j,G})}^2}},\end{array}$$(A.3)

where $g_{\mathit{MI}}^{\prime }(r_{j,G})=\frac{\partial g_{\mathit{MI}}(g_{\mathrm{bar}}(r_{j,G}))}{\partial g_{\mathrm{bar}}(r_{j,G})}$ is the derivative with respect to the baryonic acceleration.

Appendix B: Uncertainties on acceleration ratios and variance matrix

The uncertainty contributions to $\delta {\widehat{g}}_{\mathrm{bar}}(r_{j,G})$ and $\delta {\widehat{g}}_{\mathrm{obs}}(r_{j,G}),$ following from Eq. (A.1), are

$$\begin{array}{cc}\hfill \delta {\widehat{g}}_{\mathrm{obs}}(r_{j,G})& =\sqrt{{\delta }_1{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2+{\delta }_2{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2},\hfill \\ \hfill \delta {\widehat{g}}_{\mathrm{bar}}(r_{j,G})& ={\widehat{g}}_{\mathrm{bar}}(r_{j,G})\sqrt{{\left(\mathrm{\Delta }{g}_{\mathit{gas}}(r_{j,G})\right)}^2+{(\mathrm{\Delta }\mathrm{{\rm Y}}(r_{j,G}))}^2},\hfill \\ \hfill {\delta }_1{\widehat{g}}_{\mathrm{obs}}(r_{j,G})& ={\widehat{g}}_{\mathrm{obs}}(r_{j,G})\left(\frac{2\delta v_{\mathrm{obs}}(r_{j,G})}{v_{\mathrm{obs}}(r_{j,G})}\right),\hfill \\ \hfill {\delta }_2{\widehat{g}}_{\mathrm{obs}}(r_{j,G})& ={\widehat{g}}_{\mathrm{obs}}(r_{j,G})\left(\frac{2\delta v_{\mathrm{obs}}(r_{\mathrm{bar},\mathrm{G}})}{v_{\mathrm{obs}}(r_{\mathrm{bar},\mathrm{G}})}\right),\hfill \\ \hfill \mathrm{\Delta }{g}_{\mathit{gas}}(r_{j,G})& =(\frac{\delta g_{\mathrm{gas}}(r_{j,G})}{g_{\mathrm{bar}}(r_{j,G})}-\frac{\delta g_{\mathrm{gas}}(r_{\mathrm{bar},\mathrm{G}})}{g_{\mathrm{bar}}(r_{\mathrm{bar},\mathrm{G}})}),\hfill \\ \hfill \mathrm{\Delta }\mathrm{{\rm Y}}(r_{j,G})& ={\displaystyle \sum _{k=\mathrm{disc},\mathrm{bulge}}}\delta {\mathrm{{\rm Y}}}_k(\frac{v_k^2(r_{j,G})}{v_{\mathrm{bar}}^2(r_{j,G})}-\frac{v_k^2(r_{\mathrm{bar},\mathrm{G}})}{v_{\mathrm{bar}}^2(r_{\mathrm{bar},\mathrm{G}})}),\hfill \end{array}$$(B.1)

where we have separated the uncertainty contribution ${\delta }_1{\widehat{g}}_{\mathrm{obs}}(r_{j,G}),$ which is random for all points within a galaxy, and the contribution from the normalization ${\delta }_2{\widehat{g}}_{\mathrm{obs}}(r_{j,G}),$ which is a systematic for all data points within a single galaxy. With these ratios, we constructed the χ² of each galaxy compared to the model acceleration:

$$\begin{array}{c}\hfill {\chi }_G^2={\displaystyle \sum _{j,j^{\prime }}}({\widehat{g}}_{\mathrm{obs}}(r_{j,G})-{\widehat{g}}_M(r_{j,G}))V_{j{j}^{\prime }}^{-1}({\widehat{g}}_{\mathrm{obs}}(r_{j,G})-{\widehat{g}}_M(r_{j,G})),\end{array}$$(B.2)

where ${\widehat{g}}_M(r_{j,G})={\widehat{g}}_M({\widehat{g}}_{\mathrm{bar}}(r_{j,G}))$ is the MOND model prediction with the data point ${\widehat{g}}_{\mathrm{bar}}(r_{j,G})$ as input, and the inverse variance matrix (e.g. Stump et al. 2001) is

$$\begin{array}{cc}\hfill V_{j{j}^{\prime }}^{-1}& =\frac{{\delta }_{j{j}^{\prime }}}{{\delta }_1{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2}-\frac{{\delta }_2{\widehat{g}}_{\mathrm{obs}}(r_{j,G}){\delta }_2{\widehat{g}}_{\mathrm{obs}}(r_{j^{\prime },G})}{{\delta }_1{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2{\delta }_1{\widehat{g}}_{\mathrm{obs}}{(r_{j^{\prime },G})}^2}{A}^{-1},\hfill \\ \hfill A^{-1}& =1+{\displaystyle \sum _j}\frac{{\delta }_2{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2}{{\delta }_1{\widehat{g}}_{\mathrm{obs}}{(r_{j,G})}^2}.\hfill \end{array}$$(B.3)

The second term in $V_{j{j}^{\prime }}^{-1}$ comes from the systematic uncertainties ${\delta }_2{\widehat{g}}_{\mathrm{obs}}(r_{j,G})$ introduced due to the common normalization point. We neglected the uncertainties in ${\widehat{g}}_{\mathrm{bar}}(r_{\mathrm{bar},\mathrm{G}})$, which are small compared to those in ${\widehat{g}}_{\mathrm{obs}}(r_{j,G})$ for most data points, as seen in Figs. 7 and 6.

Appendix C: 3MN model

As discussed in the main text, the sum of three MN discs can be used to model an exponential disc profile using the following procedure (Smith et al. 2015). First, the single scale length used in all three potentials, b, is found in terms of the exponential disc scale length, R_d, and scale height, h_z, as

$$\begin{array}{c}\hfill \frac{b}{R_d}=-0.269{\left(\frac{h_z}{R_d}\right)}^3+1.080{\left(\frac{h_z}{R_d}\right)}^2+1.092\left(\frac{h_z}{R_d}\right).\end{array}$$(C.1)

The remaining six parameters, the three mass scale parameters, M_MN, 1, M_MN, 2, and M_MN, 3, and the three scale height parameters, a₁, a₂, and a₃, are found from the equation

$$\begin{array}{c}\hfill \text{parameter}=k_1{\left(\frac{b}{R_d}\right)}^4+k_2{\left(\frac{b}{R_d}\right)}^3+k_3{\left(\frac{b}{R_d}\right)}^2+k_4\left(\frac{b}{R_d}\right)+k_5,\end{array}$$(C.2)

with the numerical parameters k_i given in Table C.1

The 3MN models found by this table match the analytical exponential disc < 1.0% out to 4R_d and < 3.3% out to 10R_d (Smith et al. 2015). For our purposes, the 3MN model is convenient because we can use it to explicitly interpolate between spherical and disc-like matter distributions.

All Tables

All Figures

Fig. 1. Top panel: centripetal acceleration curves in g2 space with the quantities rtot, rN, and 𝒞±, used for classification in Table 1, shown. The solid grey line shows a MOND MI curve with the radius of maximum baryonic acceleration and maximum total acceleration rN = rtot indicated. Also, the curve segments 𝒞+ and 𝒞− coincide, so the curve area is 𝒜(𝒞) = 0. The dotted grey and dotted black curves show the 𝒞± curve segments of a MOND MG curve, using the Brada-Milgrom approximation in Eq. (9), with rtot >  rN and 𝒜(𝒞) > 0. For both, the baryonic matter is an infinitely thin exponential disc, Σ(r) = Σ0e−r/rd. The arrow indicates the direction of increasing radius along the curve. Bottom panel: corresponding rotation curves.

In the text

Fig. 2. Gradient of the norm of the Newtonian acceleration, ∇gN (left panel), the Newtonian acceleration vector, gN (middle panel), and the curl of the solenoidal field, ∇ × σ (right panel), in QUMOND for an MN disc model with mass parameter MMN = 1, scale height b = 0.5, and scale length a = 1. From the sign of ∇ × σ, we infer that the centripetal component of the curl acceleration is in the opposite direction of the Newtonian acceleration vector, gN, in the disc near z = 0.

In the text

Fig. 3. Top-left panel: centripetal acceleration curves in g2 space of MOND MI (solid grey curve) and MOND MG in the Brada-Milgrom approximation (black dotted curve) with acceleration scale g0 = 1.2 × 10−10 m s−2 for a 3MN model of an infinitely thin exponential disc galaxy with M = 1.2 × 1010 M⊙ and a scale length of rd = 3.5 kpc. The 3MN parameters are given in the appendix. The dotted coloured curves are the full MOND MG solution of the same model but with the scale height parameter of the 3MN model varied: b = 3 (green), b = 1 (red), and b = 0.3 (blue) (b → 0 is the disc limit, and b → ∞ is the spherical limit). Finally, the solid orange line gives the purely Newtonian acceleration. Top-right panel: corresponding rotation curves showing the full MOND MG solution only for b = 0.3. Bottom-left panel: corresponding curl field S for the b = 3 case, showing how the curl field at z = 0 is in the opposite direction to the Newtonian one and how it reduces the radial MOND MG acceleration compared to the MOND MI. Bottom-right panel: difference in curl fields for the most disc-like (b = 0.3) and most spherical (b = 3) parameter values, showing that the solenoidal field reduces the radial MOND MG acceleration in the z = 0 plane most for the most disc-like mass distributions.

In the text

Fig. 4. Top-left panel: acceleration curves of galaxies from SPARC in g2 space with cored geometry robs >  rbar. Top-right panel: same but for cuspy galaxies with robs <  rbar. Also shown in both panels are the MOND MI curves for the two interpolation functions considered in Eq. (6). Middle-left panel: corresponding model curves for the cored galaxies from MOND MG curves using the Brada-Milgrom approximation with the gbar and surface density Σ(r) values from SPARC. Middle-right panel: same as left but for the cuspy galaxies. For reference, the MOND MI curve is also shown again. Bottom-left panel: corresponding rotation speed curve data for two of the cored galaxies, compared to the MOND MI and MOND MG curves. Bottom-right panel: same as left but for two cuspy galaxies

In the text

Fig. 5. Top panel: distribution of the normalized residuals R(rj, G) in Eq. (A.3) of SPARC data with respect to MOND MI with g0 = 1.2 × 10−10 m s−2. The grey histograms with black dots at each bin top shown on both panels are the 1983 SPARC data points at large radii, i.e., rj, G >  rbar, G, after imposing the basic quality criteria and the δvobs/vobs <  0.1 cut employed in McGaugh et al. (2016), Lelli et al. (2017), and Lelli et al. (2016). Also shown in the top panel are normalized residuals of data points at small radii, rj, G ≤ rbar, G, from the (cored) galaxies only with robs, G >  rbar, G (red histogram). Bottom panel: points at small radii, rj, G <  rbar, G, from (cuspy) galaxies only with robs, G <  rbar, G (cuspy; blue histogram).

In the text

Fig. 6. Acceleration ratio curves $({\widehat{g}}_{\mathrm{bar}},{\widehat{g}}_{\mathrm{obs}})$ (points with errors) in $\widehat{g}2$ space for three of the cuspy SPARC galaxies shown in the right-hand panels of Fig. 4, compared to the MOND MI curve with g0 = 1.2 × 10−10 (solid grey), the MOND MI curve with the best-fit value of g0, min for each galaxy (dashed grey), and MOND MG in the Brada-Milgrom approximation curves (points without errors). Also shown are the χ2/d.o.f. values of the MOND MI curve with g0 = 1.2 × 10−10.

In the text

Fig. 7. $\widehat{g}2$-space data points (points with errors) of three of the cored SPARC galaxies shown in the left-hand panel of Fig. 4, compared to the MOND MI curve with g0 = 1.2 × 10−10 (solid grey), the MOND MI curve with the best-fit value of g0 for each galaxy (dashed grey), and MOND MG (Brada-Milgrom approximation) curves (points without errors). We note that most ${\widehat{g}}_{\mathrm{bar}}$ uncertainties are small for most cases because the ratios are constructed to reduce systematic uncertainties. Also shown are the χ2/d.o.f. values for the dashed best-fit MOND MI curves.

In the text