© ESO 2018

1 Introduction

Since the discovery of the flat rotation curves of disk galaxies (Bosma 1978; Rubin et al. 1978), the mass discrepancy problem has been widely explored. The baryonic Tully-Fisher relation (BTFR; Tully & Fisher 1977; McGaugh et al. 2000; Lelli et al. 2016b) was established as the link between the flat rotation velocity V _f and the baryonic mass for late-type galaxies. The definition of mass discrepancy at each radius, $M_{\text{tot}}/M_{\text{bar}}\simeq V_{\text{obs}}^2/V_{\text{bar}}^2$, makes it possible to study the “local” relation between the rotation curve shape and the baryonic mass distribution, which lead to the mass discrepancy-acceleration relation (McGaugh 2004).

In order to explore the mass discrepancy-acceleration relation further, Lelli et al. (2016a) built the Spitzer Photometry and Accurate Rotation Curves (SPARC) database: a sample of 175 disk galaxies with homogeneous [3.6] surface photometry and high-quality HI/Hα rotation curves, spanning a wide range in morphological types (S0 to Irr), stellar masses (5 dex), surface brightnesses (4 dex), and gas fractions. Using the SPARC database, McGaugh et al. (2016) established the radial acceleration relation (RAR), in which the observed acceleration ($g_{\text{obs}}=V_{\text{obs}}^2/R$) tightly correlates with the baryonic one (g_bar). The g_obs–g_bar plane has a major advantage over the mass discrepancy-acceleration relation: the two quantities and the corresponding errors are fully independent, thus observed and expected scatters can be easily computed without additional complications from covariances between the measurements. Furthermore, Lelli et al. (2017a) extend the galaxy sample to include 25 early-type galaxies and 62 dwarf spheroidals, finding that they follow the same relation as late-type galaxies within the uncertainties.

Assuming that the stellar mass-to-light ratio ϒ_⋆ does not vary strongly at [3.6] (McGaugh & Schombert 2014, 2015; Meidt et al. 2014; Schombert & McGaugh 2014), it is found that the RAR has an observed rms scatter of only 0.13 dex (McGaugh et al. 2016; Lelli et al. 2017a). This is largely driven by uncertainties on galaxy distance and disk inclination, as well as possible galaxy-to-galaxy variations in the value of ϒ_⋆. Hence, the intrinsic scatter around the RAR must be even smaller.

Given that late-type galaxies statistically satisfy the RAR, we can explore its intrinsic scatter by fitting individual galaxies and marginalizing over ϒ_⋆, galaxy distance (D), and disk inclination (i). This is equivalent to rotation curve fits in modified Newtonian dynamics (MOND; Milgrom 1983), but our aim here is to measure the intrinsic scatter around the RAR from a purely empirical perspective. Moreover, differently from classic MOND studies (e.g., Sanders & McGaugh 2002), we impose priors on ϒ_⋆, D, and i based on the observational uncertainties. These “free” parameters are treated as global quantities for each galaxy, whereas the RAR involves local quantities measured at each radius. Hence, there is no guarantee that adjusting those parameters within the errors can result in satisfactory individual fits for each and every galaxy or decrease the empirical scatter around the mean relation.

In Sect. 2, we describe our fitting method, which is a Markov chain Monte Carlo (MCMC) simulation. In Sect. 3, we show the fitted individual galaxies and their posterior distribution. The distributions of adjusted parameters are also presented. In Sect. 4, the RAR and its residuals are described. We also check the resulting BTFR. We generalize the RAR to consider possible galaxy-to-galaxy variation in the critical acceleration scale in Sect. 5. In Sect. 6, we summarize our results and discuss the general implication of the extremely small rms scatter (0.057 dex) of our best-fitting relation.

Fig. 1 Example of MCMC fits (left) and the corresponding posterior distribution (right). In the left panels, the points with error bars are the observed rotation curves V obs (R) or corresponding accelerations $V_{\text{obs}}^2(R)/R$. In the rotation curve panel, each baryonic component is presented: purple dash-dotted line for the bulge, blue dashed line for the disk, and green dotted line for the gas. The red solid line is the fitted rotation curve. The dark gray and light gray bands show the 68% and 95% confidence regions, respectively, considering the posterior distribution of ϒ⋆; they do not include additional uncertainties on i and D. In the acceleration panels, the red solid line represents the mean RAR to which we fit. In the right panels, the blue cross indicates the parameter set corresponding to the maximum posterior probability. The complete figure set (175 images) is shown in the appendix.

2 Method

2.1 Parameter dependence

Based on Spitzer [3.6] images and total HI maps from the SPARC database, one can calculate the acceleration due to the baryonic mass distribution at every radius, $$g_{\text{bar}}(R)=({{\rm Y}}_{\text{disk}}{V}_{\text{disk}}^2+{{\rm Y}}_{\text{bul}}{V}_{\text{bul}}^2+V_{\text{gas}}^2)/R,$$(1)

where ϒ_disk and ϒ _bul are the stellar mass-to-light ratios for the disk and bulge, respectively. Similarly, the observed acceleration can be calculated directly from the observed velocity V _obs, $$g_{\text{obs}}(R)=\frac{V_{\text{obs}}^2}{R}.$$(2)

According to the RAR (McGaugh 2014; McGaugh et al. 2016; Lelli et al. 2017a), the expected total acceleration g_tot strongly correlates with that expected from baryonic distributions g_bar, $$g_{\text{tot}}(R)=\frac{g_{\text{bar}}}{1-e^{-\sqrt{g_{\text{bar}}/g_{†}}}},$$(3)

where g_† = 1.20 × 10⁻¹⁰ m s⁻². Thus, one can compare the observed acceleration with the expected one.

A constant value of ϒ_⋆ for all galaxies is able to statistically establish the RAR, but some scatter must have been introduced since ϒ _⋆ should vary from galaxy to galaxy. An inappropriate ϒ_⋆ can lead to systematic offsets from the RAR for individual objects. Specifically, ϒ _disk and ϒ _bulge affect g_bar according toEq. (1).

Uncertainties in galaxy distance affect the radius (R) and the baryonic components of the rotation curve (V _k). With D being adjusted to D′, R and V _k transform as $$R^{\prime }=R\frac{D^{\prime }}{D};{V^{\prime }}_k=V_k\sqrt{\frac{D^{\prime }}{D}},$$(4)

where k denotes disk, bulge or gas. Therefore, g_bar does not depend on distance. Instead, g_obs goes as D⁻¹ because the observed rotation velocity (V _obs) and its error (δV _obs) are inferred from the line-of-sight velocity which is distance independent.

In the SPARC database, galaxy distances are estimated using five different methods (see Lelli et al. 2016a for details): (1) the Hubble flow corrected for Virgo-centric infall (97 galaxies), (2) the tip of the red giant branch (TRGB) method (45 galaxies), (3) the magnitude-period relation of Cepheids (3 galaxies), (4) membership to the Ursa Major cluster of galaxies (28 galaxies), and (5) supernovae (SN) light curves (2 galaxies). The first method is the least accurate because the systemic velocity (redshift) of a galaxy may be largely affected by peculiar flows in the nearby Universe. The other methods have accuracies ranging between 5% and 15%. Table 1 in McGaugh et al. (2016) shows that errors on galaxy distance are the main source of scatter on the RAR.

Uncertainties in disk inclination are another important source of scatter. When the disk inclination i is adjusted to i′, V _obs and δV _obs transform as $${V^{\prime }}_{\text{obs}}=V_{\text{obs}}\frac{\sin (i)}{\sin (i^{\prime })};\delta {V^{\prime }}_{\text{obs}}=\delta V_{\text{obs}}\frac{\sin (i)}{\sin (i^{\prime })}.$$(5)

Hence, g_obs has a further dependence on disk inclination. Clearly, the correction becomes very large for face-on galaxies with small inclination. We also note that several galaxies have warped HI disks: the inclination angle systematically varies with radius. While warps are taken into account in deriving rotation curves (e.g., NGC 5055 from Battaglia et al. 2006), here we treat the inclination as a single globalparameter for each galaxy.

2.2 MCMC simulation

To fit individual galaxies, we used emcee (Foreman-Mackey et al. 2013) to map the posterior distribution of the parameter set: the stellar mass-to-light ratio, galaxy distance, and disk inclination. Following standard procedures, we imposed Gaussian priors on ϒ _⋆, D, and i. The priors were centered around the assumed values in SPARC and have standard deviations given by the observational errors for D and i and the scatter expected from stellar population models for ϒ_⋆ (e.g., Bell & de Jong 2001). Hence, ϒ_⋆, D, and i are not entirely free parameters: the MCMC simulation is searching for an optimal solution within a realistic region of the parameter space. Specifically, we imposed ϒ_disk = 0.5 and ϒ_bulge = 0.7 M_⊙∕L_⊙ with a standard deviation of 0.1 dex. We adopted a fixed mass-to-light conversion for the gas unlike what Swaters et al. (2012) did. We also required that the parameters remain physical and positive definite: ϒ _⋆ > 0, D > 0 Mpc, and 0° < i < 90°.

We used the standard affine-invariant ensemble sampler in emcee and initialized the MCMC chains with 200 random walkers. We ran 500 burnt-in iterations and then ran the simulation to more than five autocorrelation times. We checked that the acceptance fractions for all galaxies are in the range (0.1, 0.7). To achieve the acceptance fraction, we set the size of the stretch move a = 2.

We record the parameter set corresponding to the maximum probability and calculate the reduced χ², $${\chi }_{\nu }^2={\sum }_R\frac{{[g_{\text{obs}}(R)-g_{\text{tot}}(R)]}^2/{\sigma }_{g_{\text{obs}}}^2}{N-f},$$(6)

where ${\sigma }_{g_{\text{obs}}}=2V_{\text{obs}}\times \frac{\delta V_{\text{obs}}}{R}$ is the uncertainty in the observedacceleration, N the number of data points, and f the degrees of freedom, for every galaxy.

Fig. 2 Same as Fig. 1 but for the gas-dominated dwarf galaxy IC 2574.

3 Fittingindividual galaxies

3.1 Fit results

By fitting individual galaxies to the RAR, the stellar mass-to-light ratio, the galaxy distance and the disk inclination are optimized. In particular, we note that the RAR may be used as a distance indicator in analogy to the BTFR. The former relation is more demanding in terms of data quality, but has the advantage of using the full shape of the rotation curve and the baryonic mass profile instead of merely using the flat rotation velocity and total baryonic mass that go into the BTFR.

Figure 1 shows an example of an MCMC fit for a star-dominated spiral galaxy (NGC 2841). This object has historically been regarded as a problematic case for MOND (Begeman et al. 1991; Gentile et al. 2011), but a good fit is obtained allowing for uncertainties in distance and inclination within 1σ. The values of ϒ_disk and ϒ_bulge are relatively high but not unreasonable for such a massive, metal-rich galaxy. Similar figures are available for all SPARC galaxies.

Figure 2 illustrates the MCMC fit of a gas-dominated dwarf galaxy (IC 2574). This object is often considered a problematic case for Λ cold dark matter (CDM) because it has a large core extending over ~8 kpc (e.g., Oman et al. 2015). Moreover, Navarro et al. (2017) claimed that this galaxy strongly deviates from the RAR (their Fig. 3). We find an excellent fit for IC 2574 after adjusting its distance and inclination by 1σ_D and 1.5σ_i, respectively. The adjusted mass-to-light ratio is rather low (ϒ_disk = 0.07 M_⊙ ∕L_⊙), perhaps uncomfortably so. This object is also present in the THINGS database (de Blok et al. 2008), where the rotation curve is consistent with but slightly higher than that adopted here. If we apply the same MCMC technique to the THINGS data, we find a good fit with ϒ_⋆ = 0.25 M_⊙ ∕L_⊙, illustrating the sensitivity of this parameter to even small changes in the rotation curve.

We stress that for gas-rich dwarfs, the vast majority of the rotation curve is explained by the gas contribution with very little room for adjustment. Rather than be overly concerned with the exact value of ϒ _⋆ in such gas-dominated galaxies, it is amazing that this procedure works at all: ϒ _⋆ has little power to affect the overall fit, while D and i are constrained by their priors to be consistent with the observed values. Gas-dominated galaxies are more prediction than fit: given the observed gas distribution and the RAR, the rotation curve must be what it is. The fitting parameters provide only minor tweaks to the basic prediction.

In general, the fits to most galaxies are good. The mass-to-light ratios are generally consistent with the expectations of stellar population synthesis. It is rare that either D or i are adjusted outside of their observational uncertainties. We maintain the same fitting function (Eq. (3)) for all 175 galaxies.

While most fits are visually good, they may occasionally have poor values of χ². These usually occur when one or a few individual velocity measurements have tiny error bars. The discrepancy of these points from the fit is small in an absolute sense, but still impacts χ². It is likely that in some cases the errors are slightly underestimated. For example, the potential contribution of non-circular motions may have been understated and the velocities may not exactly trace the underlying gravitational potential. In general, these fits are as good as possible: one cannot do better with a dark matter halo fit. The Navarro-Frenk-White (NFW) halo fit to NGC 2841 (Katz et al. 2017) looks indistinguishable from that in Fig. 1: the two extra fit parameters available with a dark matter fit do not alter the shape of the continuous line that best approximates the data. We therefore consider fits of this type to be good even if χ² is larger than unity.

In about 10% of the cases, however, the fits are genuinely poor. Poor fits generally happen for rotation curves of lowest quality. For the sake of completeness, we fit all 175 galaxies in the SPARC database, but we did not expect to find good fits for galaxies with quality flag Q = 3 (e.g., NGC 4389, UGC 2455, UGC 4305) where the gas kinematics is likely out of equilibrium (see Lelli et al. 2016b, for details). Some poor fits are also found for galaxies with Q = 1 (e.g., D631-7, F571-8, and IC 4202) and Q = 2 (e.g., Cam B, DDO 168, and NGC 2915). This may happen for several reasons: (i) the errors on D and i may be slightly underestimated, hence the priors place strong constraints that then contribute more to the probability function and preclude somewhat better fits, (ii) there may be features in the rotation curves that do not tracethe smooth gravitational potential but are due to large-scale non-circular motions, and (iii) these galaxies may have unusual dust content that affects the shape of the [3.6] luminosity profile and the calculation of g_bar (this may be particularly important for edge-on systems such as IC 4202 and F571-8). In the specific case of D631-7, the gas contribution was computed assuming a purely exponential distribution since the HI surface density profile was not available (see Lelli et al. 2016a). Since the RAR is very sensitive to the precise baryonic distribution, even small deviations from an exponential profile may lead to a poor fit in such a gas-dominated galaxy. In general, we consider it likely that lower quality data lead to lower quality fits. Having ~10% of such cases seems an inevitable occurrence in any astronomical database built from many diverse rotation curve studies such as SPARC (Lelli et al. 2016a). It would be strange if there were no such cases.

Fig. 3 Distributions of fitted parameters. The top panels show the histograms of the optimized ϒ disk (top left) and ϒbulge (top right). The vertical dashed lines represent the values of 0.5 (disk) and 0.7 (bulge) adopted in McGaugh et al. (2016). In the bottom panels, the optimal galaxy distance (bottom left) and disk inclination (bottom right) are plotted against their original values. The dashed line is the line of unity. Different methods of measuring distance are indicated by different colors. Large and small symbols correspond to data with an accuracy higher and lower than 15% for distance and 5% for inclination, respectively, based on observational errors as tabulated in SPARC. Crosses indicate galaxies with the low-quality rotation curves (Q = 3, see Lelli et al. 2016a). A few other outliers are discussed in the text.

Fig. 4 RAR (left) and the residuals (right) with ϒ ⋆, D, and i optimized by the MCMC method. In the left panel, the red solid line represents the mean RAR from Eq. (3). The black dotted line is the line of unity. 2694 individual data points from 153 SPARC galaxies are represented by the blue color-scale. Black dashed lines show the rms scatter. In the right panel, the histogram of the residuals is fit with both single (black solid line) and double (red solid line) Gaussian functions. Blue dashed lines show the two components of the double Gaussian function.

3.2 Distributions of adjusted parameters

Figure 3 shows the distributions of the optimized parameters. The top panels show histograms of ϒ _⋆. The dashed lines indicate ϒ_disk = 0.5 M_⊙ ∕L_⊙ and ϒ _bulge = 0.7 M_⊙ ∕L_⊙ adopted in McGaugh et al. (2016) and Lelli et al. (2017a). The optimized stellar mass-to-light ratios are tightly distributed around these values. The median values of ϒ_disk and ϒ_bulge are 0.50 and 0.67, respectively. By and large, the best-fit stellar mass-to-light ratios are consistent with the expectations of stellar population synthesis models (Schombert & McGaugh 2014; Meidt et al. 2014; McGaugh & Schombert 2014; Norris et al. 2016).

Adjusted distances and inclinations are also shown in Fig. 3 (bottom panels). Galaxies with a low-quality flag in SPARC (Q = 3, see Lelli et al. 2016a) typically prefer smaller values of D and i with respect to their original values (see black crosses in the figure). After removing these low-quality data, the distributions of D and i are fairly symmetric around the line of unity indicating that there are no major systematics. Hubble flow distances are the least certain and show the largest variation between measured and best-fit distance. More accurate methods (Cepheids, TRGB) show less variation, as expected.

A few galaxies show significant deviations (>1σ_D) from the optimized ones. For example, PGC51017 (with a TRGB distance) is a starburst dwarf galaxy where the rotation curve clearly does not trace the equilibrium gravitational potential (Lelli et al. 2014), thus it is not surprising that the distance is pushed to unphysical values in order to obtain a good fit. This rotation curve has a low-quality flag in SPARC (Q = 3) and is only included here for the sake of completeness. Another example is NGC3198, which is sometimes regarded as a problematic case for MOND (Gentile et al. 2011, 2013). The MCMC method finds a good fit with D = 10.4 ± 0.4 Mpc, which is consistent with the Cepheid-based distance (13.8 ± 1.4 Mpc) within 2σ.

Table A.1 lists the optimal parameters and reduced χ² for each galaxy in order of declining luminosity. The errors on the fitting parameters are estimated from their posterior distributions using the “std” output in GetDist. These errors are generally smaller than those in the SPARC database because of the combined constraints from the Gaussian priors and the likelihood function.

4 Galaxy scaling relations

4.1 Radial acceleration relation

The RAR and its residuals, log(g_obs) − log(g_tot), are plotted in Fig. 4. To compare them with previous results (McGaugh et al. 2016; Lelli et al. 2017a), the same selection criteria were adopted: we removed 10 face-on galaxies with i < 30° and 12 galaxies with asymmetric rotation curves that do not trace the equilibrium gravitational potential (Q = 3). We also required a minimum precision of 10% in observational velocity (δV _obs∕V _obs < 0.1). This retains 2694 data points out of 3163.

After ϒ_⋆, D, and i were adjusted within the errors, the RAR was extremely tight and had an rms scatter of 0.057 dex. We fit the histogram of the residuals with a Gaussian function (the dashed line in the right panel of Fig. 4): the fit is acceptable, but there are broad symmetric wings in the residuals that are not captured by a single Gaussian function. Hence, we fit a double Gaussian function. The double Gaussian function substantially improves the fit and fully describes the residual distribution with standard deviations of 0.062 dex and 0.020 dex. The mean values μ are consistent with zero.

Interestingly, the errors on the rotation velocities are not expected to be Gaussian because they are obtained by summing two different contributors (see Swaters et al. 2009; Lelli et al. 2016a): the formal error from fitting the whole disk (driven by data quality and random non-circular motions) and the difference between velocities in the approaching and receding sides of the galaxy (representing global asymmetries and kinematic lopsidedness). The success of the double Gaussian fit suggests that the two Gaussian components perhaps probe these two different sources of errors and hence dominate the total residual scatter over all other possible error sources.

The small residual scatter leaves very little room for any intrinsic scatter because (1) the observational errors in the rotation velocities are not negligible, driving errors in g_obs, and (2) there could be additional sources of errors in g_bar like the detailed 3D geometry of baryons and possible radial variations in ϒ _⋆. Considering these error sources, the intrinsic scatter in the RAR must be smaller than 0.057 dex.

4.2 Baryonic Tully-Fisher relation

The BTFR relates the total baryonic masses of galaxies to their flat rotation velocity. In some sense, this is the asymptotic version of the RAR at large radii (see Sect. 7.1 of Lelli et al. 2017a for details). For R →∞, g_bar becomes small and Eq. (3) gives $g_{\text{obs}}\simeq \sqrt{g_{\text{bar}}{g}_{†}}$ by Taylor expansion. Since $g_{\text{obs}}=V_{\text{f}}^2/R$ and g_bar ≃ GM_bar∕R², the radial dependence cancels out and we are left with $V_{\text{f}}^4\propto M_{\text{bar}}$. Thus, a BTFR with slope 4 is built into Eq. (3).

Here, we fit the BTFR directly to check how well we recover the behavior required by Eq. (3). In addition to the slope specified by the asymptotic limit of the RAR, we also expect the BTFR to be tighter after fitting ϒ _⋆, D, and i to the RAR. However, this does not necessarily have to happen since the BTFR only considers the flat rotation velocity (V _f) and the total baryonic mass (M_b), whereas we fit the whole shape of the rotation curve using the full baryonic mass profile. We adopted the same selection criteria as described in Sect. 4.1 and removed the galaxies that did not reach a flat rotation velocity as defined by Lelli et al. (2016b). This retained 123 galaxies out of 175 (5 more galaxies with the latest version of SPARC relative to Lelli et al. 2016b).

Figure 5 shows a tight BTFR. We used the LTS_LINEFIT program (Cappellari et al. 2013) to fit the linear relation $$\log (M_{\text{b}})=n\log (V_{\text{f}})+\log (A).$$(7)

LTS_LINEFIT considers errors in both variables and estimates n and A together with the intrinsic scatter around the linear relation. The errors on V _f and M_b were calculated using Eqs. (3) and (5) of Lelli et al. (2016b), but we treated disk and bulge separately. These equations consider the errors on ϒ_disk, ϒ _bulge, i, and D, which were estimated from the posterior distributions of the MCMC fits. We also corrected observed quantities such as luminosity and flat rotation velocity according to the adjusted D and i. The fitting results are summarized in Table 1. To enable a direct comparison with Lelli et al. (2016b), we also show the case where D and i were kept fixed to the SPARC values and the mass-to-light ratio was constant for all galaxies, but we improved compared to Lelli et al. (2016b) by considering disk (ϒ_⋆ = 0.5) and bulge (ϒ _⋆ = 0.7) separately. In this case, the errors on D and i were taken from the SPARC database, while the errors on ϒ_disk and ϒ_bulge were assumed to be 0.11 dex for all galaxies. We find a slightly steeper slope than Lelli et al. (2016b) because bulges have higher mass-to-light ratios than disks and are more common in massive galaxies, increasing M_b at the top end of the relation. Except for this small difference, our results are entirely consistent with those of Lelli et al. (2016b) for the constant ϒ _⋆ case.

The estimated intrinsic scatter of the BTFR is rather small. We estimate a conservative upper limit on the intrinsic scatter as the best-fit value (0.035 dex) plus 3σ (the error on the estimated intrinsic scatter). This gives σ_intr < 0.1 dex. Satisfying this bound provides a strong constraint on galaxy formation models (e.g., Desmond 2017a).

The fitted slope of the BTFR is close to 3.8 in both cases. Formally, this differs from 4.0 by ~ 4σ, although slopes consistent with 4 are obtained when we weight the data by the gas fraction (Lelli et al. 2016b). Given that the functional form of Eq. (3) guarantees a BTFR with a slope of 4, this discrepancy is puzzling.

Several effects may cause the difference. One is simply that there are uncertainties in the mass-to-light ratios estimated from our fits. This adds scatter to the data, which inevitably lowers the fitted slope.

A more subtle concern is that the measured value of V _f is not quitethe same as that implicit in Eq. (3). We measure V _f in the outer parts of extended rotation curves, and can do so consistently and robustly. However, Eq. (3) only guarantees a BTFR slope of 4 with V _f in the limit of zero acceleration, or infinite distance from the galaxy. The measurements are made at finite radii. The definition of V _f in (Lelli et al. 2016b, their Eq. (2)) requires measured rotation curves to be flat to within 5%, but there may be some small slope within that limit. It is well known (e.g., Verheijen 2001) that bright galaxies have rotation curves that tend to decline toward V _f, while those of faint galaxies tend to rise toward V _f. It is conceivable that this effect causes a slight systematic variation in the measured V _f with mass that acts to lower the slope. That is to say, the value of V _f we measure empirically may not reach the flat velocity sufficiently well that is implied by the limit g_bar→ 0 assumed in the derivation above.

The geometry of disk galaxies may also have an impact: a thin disk rotates faster than the equivalent spherical distribution (Binney & Tremaine 1987). This is quantified by the factor X that appears in the normalization of the MOND prediction for the BTFR: A = X∕(a₀G). The factor X → 1 as R →∞, but on average, ⟨X⟩ = 0.8 (McGaugh 2005) at the finite radii observed in spiral galaxies. We have assumed that X is the same for all galaxies, but it is conceivable that disk thickness varies with mass so that X is a weak function thereof. This might also affect the slope.

Regardless of which of these effects dominates, it is clear that the slope of the BTFR is steep. It is not 3.0 as one might reasonably assume in Λ CDM (e.g., Mo et al. 1998), nor is it 3.5 (e.g., Bell & de Jong 2001), as might be expected after adiabatic compression (Bullock et al. 2001). The difference (or lack thereof) between 3.8 and 4.0 may be the limit of what we can hope to discern with astronomical data. The limit is not due to the data themselves, but to systematic effects.

Fig. 5 Baryonic Tully-Fisher relation with ϒ⋆, D, and i optimized (red points). The solid line illustrates the fitted BTFR. For reference, the dashed line shows the prediction of MOND: log(Mb) = 4log(V f) + log[X∕(a0G)].

Fig. 6 Cumulative distributions of reduced χ2 for fixed g† fits (red line)and variable g† fits using flat (black line) and Gaussian priors (blue line).

5 Does the critical acceleration scale g_† vary?

In the previous analysis we have assumed that g_† is constant for all galaxies. The small scatter observed around the RAR (McGaugh et al. 2016; Lelli et al. 2017a) already demonstrates that this is very nearly the case. However, the answer to the question of whether the value of g_† is truly constant can be used to distinguish between a scaling relation and a law of nature.

As a further check on this point, we fit all galaxies again, treating g_† as an additional free parameter. We made fits with both a flat prior (for the range 0 ≤ g_†≤ 10⁻⁹ m s⁻²) and a Gaussian prior (with g_† = 1.20 ± 0.02 m s⁻²: McGaugh et al. 2016; Lelli et al. 2017a). The cumulative distributions of reduced χ² of these fits are shown along with that for fixed g_† in Fig. 6.

Allowing g_† to vary from galaxy to galaxy does not improve the fits. In the case of the flat prior, the cumulative distribution of reduced χ² is practically indistinguishable from the case of fixed g_†, despite the additional freedom. In the case of a Gaussian prior, the reduced χ² is even slightly worse because essentially the same fit is recovered (for a similar total χ²), but the extra parameter increases the number of degrees of freedom, increasing the reduced χ². The fits are not meaningfully improved byallowing g_† to vary.

The resulting rms scatter remains nearly invariant: 0.054 dex and 0.057 dex when using the flat prior and the Gaussian prior on g_†, respectively.This indicates that the remaining rms scatter is dominated by observational uncertainties on rotation curves and possibleintrinsic scatter. As a practical matter, there is no room to accommodate substantial variation in g_†.

We show the distributions of best-fit g_† in Fig. 7 for both flat and Gaussian priors. The flat prior leads to a wide distribution of g_†, while a Gaussian prior results in a tight distribution around its fiducial value. The Gaussian prior indeed results in a distribution so close to a fixed g_†, with a width smaller than the standard deviation imposed by the prior, that it appears as a δ-function on a scale that accommodates the wide distribution of the flat prior. That the apparent distribution of g_† is so large in the case of the flat prior is indicative of parameter degeneracy: changes in g_† can be compensated for by changes in the mass-to-light ratio (or the nuisance parameters) so that both may vary in an unphysical way to achieve trivial gains in χ².

Adjusting the value of g_† improves neither the fits nor the rms scatter. The data are consistent with the same value for all galaxies. There is no need to invoke variable g_† (cf. Bottema & Pestaña 2015). To do so would violate the law of parsimony (Occam’s razor).

Fig. 7 Distribution of optimal g† imposing a flat (light blue) and Gaussian (dark blue) prior. The red dashed line marks our fiducial value 1.2 × 10−10 m s−2. The inset shows the Gaussian prior alone, zoomed in and switched to a linear scale to resolve the distribution. Note the vast difference in scales: the flat prior results in a broad range of g† (although with no improvement in χ2, as seen in Fig. 6), while the Gaussian prior effectively returns a constant g†.

6 Discussion and conclusion

By fitting individual galaxies with the MCMC method, we showed that the intrinsic scatter of the RAR is extremely small. The baryonic matter distribution can reproduce the rotation curve very well, and vice versa.

The tightness of the RAR provides a challenge for the standard Λ CDM cosmology. Recent studies claim that the RAR is a natural product of galaxy formation in Λ CDM (Keller & Wadsley 2017; Ludlow et al. 2017), but none of these studies have properly taken into account observational effects when comparing theory and observations. There are two major issues: (1) general confusion between the concepts of observed and intrinsic scatter, and (2) oversampling of simulated rotation curves. Keller & Wadsley (2017) analyzed 32 galaxies from the MUGS2 “zoom-in” hydrodynamic simulations and argued that the dissipative collapse of baryons can result in a relation with a scatter of 0.06 dex. Similarly, Ludlow et al. (2017) analyzed a suite of simulated galaxies from the EAGLE project and found a relation with a scatter of 0.09 dex. However, these values cannot be compared with the observed scatter because (1) measurement errors are not added to the simulated galaxies and are not properly propagated, and (2) the simulated rotation curves are not sampled with the same number of resolution elements as in the observations. Both effects can significantly underpredict the scatter expected from cosmological simulations. Oversampling the same error-free simulated galaxy over and over can artificially decrease the expected scatter around the mean relation.

In contrast, Desmond (2017b) took both radial sampling and observational errors into account (see also Di Cintio & Lelli 2016). Desmond (2017b) found that fiducial ΛCDM models significantly overpredict the observed scatter around the mass discrepancy-acceleration relation by ~3.5σ. This discrepancy remains even if one assumes a perfect 1:1 relation between halo mass and stellar mass of galaxies with no scatter. Hence, this problem seems to be due to the stochastic hierarchical formation of DM halos, independent of baryonic physics.

The detailed shape of the RAR is also important. Ludlow et al. (2017) fit the Eq. (3) to their simulated galaxies, but found a value of g_† = 2.6 × 10⁻¹⁰ m s⁻² instead of g_† = 1.2 ± 0.02 × 10⁻¹⁰ m s⁻². This is a 70σ discrepancy. Recently, Tenneti et al. (2018) analyzed galaxies from the MassiveBlack-II simulation. They also found a correlation between g_tot and g_bar, but this is better fit by a power law with a width of ~0.1 dex rather than by Eq. (3). Hence, the detailed properties of the RAR (shape and scatter) remain an open issue for Λ CDM models of galaxy formation.

We here reported an rms scatter of 0.057 dex after marginalizing over the uncertainties due to mass-to-light ratio, galaxy distance, and disk inclination. This observed scatter is a hard upper limit to the intrinsic scatter since the errors on the observed velocities are non-negligible. This is hard to understand in a Λ CDM scenario since the diverse formation histories of galaxies must necessarily introduce some scatter on the relation, as demonstrated by the existing cosmological simulations.

In order to reconcile the conundrum that the standard CDM faces, some new dark sectors have been proposed, such as dipolar DM particles subjected to gravitational polarization (Blanchet & Le Tiec 2009, 2008), dark fluids (Zhao & Li 2010; Khoury 2015), dissipative DM particles (Chashchina et al. 2017), or fifth forces (Burrage et al. 2017). To be viable, such hypotheses must explain the shape, amplitude, and negligible scatter of the RAR. The coupling between the baryonic matter and dark matter must be rather strong to explain these observations.

On the other hand, the tight RAR could be easily understood in MOND (Milgrom 1983). MOND dictates that the equations of motions become scale-invariant at accelerations a < a₀, where a₀ corresponds to g_† in Eq. (3) (Milgrom 2009). Thus, Eq. (3) is related to the interpolation function ν(g_bar∕a₀) of MOND. The scale invariance can be achieved in two ways: modified gravity (MG) by changing the Poisson equation (Bekenstein & Milgrom 1984; Milgrom 2010), and modified inertia (MI) by changing the Second Law of Newton (Milgrom 1994). MI requires the relation g_obs = ν(g_bar∕a₀)g_bar to be true for circular orbits, leading to zero intrinsic scatter in the RAR of rotating disk galaxies (Milgrom 1994). MG requires the system to be spherically symmetric to hold precisely to the acceleration relation (Bekenstein & Milgrom 1984). Thus, the predicted g_tot of disk galaxies can show subtle differences in MG (Brada & Milgrom 1995; Milgrom 2012), and some non-zero intrinsic scatter could be introduced. Our results suggest a preference for the MI theory.

Other modified gravity theories such as emergent gravity could also potentially explain the RAR (Verlinde 2017). However, Lelli et al. (2017b) shows that the RAR predicted by this theory has significant intrinsic scatter and the residuals should correlate with radius, which contradicts the data.

The extremely small intrinsic scatter of the RAR provides a tool for testing various theories of modified gravity or dark matter. It also provides key insights toward the path for finally solving the “dark matter problem”.

Acknowledgements

We thank Harley Katz for useful discussion about using MCMC and GetDist. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This publication was made possible in part through the support of the John Templeton Foundation.

Appendix A Additional material

Fig. A.1 Complete figure set for 175 SPARC galaxies.

References

All Tables

All Figures

Fig. 1 Example of MCMC fits (left) and the corresponding posterior distribution (right). In the left panels, the points with error bars are the observed rotation curves V obs (R) or corresponding accelerations $V_{\text{obs}}^2(R)/R$. In the rotation curve panel, each baryonic component is presented: purple dash-dotted line for the bulge, blue dashed line for the disk, and green dotted line for the gas. The red solid line is the fitted rotation curve. The dark gray and light gray bands show the 68% and 95% confidence regions, respectively, considering the posterior distribution of ϒ⋆; they do not include additional uncertainties on i and D. In the acceleration panels, the red solid line represents the mean RAR to which we fit. In the right panels, the blue cross indicates the parameter set corresponding to the maximum posterior probability. The complete figure set (175 images) is shown in the appendix.

In the text

	Fig. 2 Same as Fig. 1 but for the gas-dominated dwarf galaxy IC 2574.
In the text

Fig. 3 Distributions of fitted parameters. The top panels show the histograms of the optimized ϒ disk (top left) and ϒbulge (top right). The vertical dashed lines represent the values of 0.5 (disk) and 0.7 (bulge) adopted in McGaugh et al. (2016). In the bottom panels, the optimal galaxy distance (bottom left) and disk inclination (bottom right) are plotted against their original values. The dashed line is the line of unity. Different methods of measuring distance are indicated by different colors. Large and small symbols correspond to data with an accuracy higher and lower than 15% for distance and 5% for inclination, respectively, based on observational errors as tabulated in SPARC. Crosses indicate galaxies with the low-quality rotation curves (Q = 3, see Lelli et al. 2016a). A few other outliers are discussed in the text.

In the text

Fig. 4 RAR (left) and the residuals (right) with ϒ ⋆, D, and i optimized by the MCMC method. In the left panel, the red solid line represents the mean RAR from Eq. (3). The black dotted line is the line of unity. 2694 individual data points from 153 SPARC galaxies are represented by the blue color-scale. Black dashed lines show the rms scatter. In the right panel, the histogram of the residuals is fit with both single (black solid line) and double (red solid line) Gaussian functions. Blue dashed lines show the two components of the double Gaussian function.

In the text

	Fig. 5 Baryonic Tully-Fisher relation with ϒ⋆, D, and i optimized (red points). The solid line illustrates the fitted BTFR. For reference, the dashed line shows the prediction of MOND: log(Mb) = 4log(V f) + log[X∕(a0G)].
In the text

	Fig. 6 Cumulative distributions of reduced χ2 for fixed g† fits (red line)and variable g† fits using flat (black line) and Gaussian priors (blue line).
In the text

Fig. 7 Distribution of optimal g† imposing a flat (light blue) and Gaussian (dark blue) prior. The red dashed line marks our fiducial value 1.2 × 10−10 m s−2. The inset shows the Gaussian prior alone, zoomed in and switched to a linear scale to resolve the distribution. Note the vast difference in scales: the flat prior results in a broad range of g† (although with no improvement in χ2, as seen in Fig. 6), while the Gaussian prior effectively returns a constant g†.

In the text

	Fig. A.1 Complete figure set for 175 SPARC galaxies.
In the text