Issue 
A&A
Volume 615, July 2018



Article Number  A3  
Number of page(s)  70  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201732547  
Published online  03 July 2018 
Fitting the radial acceleration relation to individual SPARC galaxies
^{1}
Deparment of Astronomy, Case Western Reserve University,
Cleveland
OH 44106, USA
email: pengfeili0606@gmail.com
^{2}
European Southern Observatory,
KarlSchwarschildStrasse 2,
Garching bei München, Germany
^{3}
Department of Physics, University of Oregon,
Eugene
OR 97403, USA
Received:
24
December
2017
Accepted:
28
February
2018
Galaxies follow a tight radial acceleration relation (RAR): the acceleration observed at every radius correlates with that expected from the distribution of baryons. We use the Markov chain Monte Carlo method to fit the mean RAR to 175 individual galaxies in the SPARC database, marginalizing over stellar masstolight ratio (ϒ_{⋆}), galaxy distance, and disk inclination. Acceptable fits with astrophysically reasonable parameters are found for the vast majority of galaxies. The residuals around these fits have an rms scatter of only 0.057 dex (~13%). This is in agreement with the predictions of modified Newtonian dynamics (MOND). We further consider a generalized version of the RAR that, unlike MOND, permits galaxytogalaxy variation in the critical acceleration scale. The fits are not improved with this additional freedom: there is no credible indication of variation in the critical acceleration scale. The data are consistent with the action of a single effective force law. The apparent universality of the acceleration scale and the small residual scatter are key to understanding galaxies.
Key words: dark matter / Galaxy: kinematics and dynamics / galaxies: spiral / galaxies: dwarf / galaxies: irregular
© 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 TullyFisher 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 latetype 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 discrepancyacceleration relation (McGaugh 2004).
In order to explore the mass discrepancyacceleration 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 highquality 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 discrepancyacceleration 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 earlytype galaxies and 62 dwarf spheroidals, finding that they follow the same relation as latetype galaxies within the uncertainties.
Assuming that the stellar masstolight 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 galaxytogalaxy variations in the value of ϒ_{⋆}. Hence, the intrinsic scatter around the RAR must be even smaller.
Given that latetype 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 galaxytogalaxy 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 bestfitting 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 dashdotted 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 masstolight 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}_{\u2020}}}},$$(3)
where g_{†} = 1.20 × 10^{−10} m s^{−2}. 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};\text{\hspace{0.17em}\hspace{0.17em}}{{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^{−1} because the observed rotation velocity (V _{obs}) and its error (δV _{obs}) are inferred from the lineofsight 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 Virgocentric infall (97 galaxies), (2) the tip of the red giant branch (TRGB) method (45 galaxies), (3) the magnitudeperiod 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{\mathrm{sin}(i)}{\mathrm{sin}({i}^{\prime})};\text{\hspace{0.17em}\hspace{0.17em}}\delta {{V}^{\prime}}_{\text{obs}}=\delta {V}_{\text{obs}}\frac{\mathrm{sin}(i)}{\mathrm{sin}({i}^{\prime})}.$$(5)
Hence, g_{obs} has a further dependence on disk inclination. Clearly, the correction becomes very large for faceon 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 (ForemanMackey et al. 2013) to map the posterior distribution of the parameter set: the stellar masstolight 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 masstolight 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 affineinvariant ensemble sampler in emcee and initialized the MCMC chains with 200 random walkers. We ran 500 burntin 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 χ^{2}, $${\chi}_{\nu}^{2}={\sum}_{R}\frac{{[{g}_{\text{obs}}(R){g}_{\text{tot}}(R)]}^{2}/{\sigma}_{{g}_{\text{obs}}}^{2}}{Nf},$$(6)
where ${\sigma}_{{g}_{\text{obs}}}=2{V}_{\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.
3 Fittingindividual galaxies
3.1 Fit results
By fitting individual galaxies to the RAR, the stellar masstolight 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 stardominated 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, metalrich galaxy. Similar figures are available for all SPARC galaxies.
Figure 2 illustrates the MCMC fit of a gasdominated 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 masstolight 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 gasrich 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 gasdominated 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. Gasdominated 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 masstolight 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 χ^{2}. 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 χ^{2}. It is likely that in some cases the errors are slightly underestimated. For example, the potential contribution of noncircular 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 NavarroFrenkWhite (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 χ^{2} 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., D6317, F5718, 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 largescale noncircular 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 edgeon systems such as IC 4202 and F5718). In the specific case of D6317, 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 gasdominated 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 lowquality 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 colorscale. 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 masstolight 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 bestfit stellar masstolight 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 lowquality 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 lowquality 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 bestfit 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 lowquality 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 Cepheidbased distance (13.8 ± 1.4 Mpc) within 2σ.
Table A.1 lists the optimal parameters and reduced χ^{2} 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 faceon 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 noncircular 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.
BTFR: the fitted parameters and scatter.
4.2 Baryonic TullyFisher 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}_{\u2020}}$ by Taylor expansion. Since ${g}_{\text{obs}}={V}_{\text{f}}^{2}/R$ and g_{bar} ≃ GM_{bar}∕R^{2}, 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 $$\mathrm{log}({M}_{\text{b}})=n\mathrm{log}({V}_{\text{f}})+\mathrm{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 masstolight 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 masstolight 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 bestfit 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 masstolight 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_{0}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 TullyFisher 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(M_{b}) = 4log(V _{f}) + log[X∕(a_{0}G)]. 
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^{−9} m s^{−2}) and a Gaussian prior (with g_{†} = 1.20 ± 0.02 m s^{−2}: McGaugh et al. 2016; Lelli et al. 2017a). The cumulative distributions of reduced χ^{2} 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 χ^{2} is practically indistinguishable from the case of fixed g_{†}, despite the additional freedom. In the case of a Gaussian prior, the reduced χ^{2} is even slightly worse because essentially the same fit is recovered (for a similar total χ^{2}), but the extra parameter increases the number of degrees of freedom, increasing the reduced χ^{2}. 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 bestfit 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 masstolight ratio (or the nuisance parameters) so that both may vary in an unphysical way to achieve trivial gains in χ^{2}.
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 “zoomin” 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 errorfree 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 discrepancyacceleration 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^{−10} m s^{−2} instead of g_{†} = 1.2 ± 0.02 × 10^{−10} m s^{−2}. This is a 70σ discrepancy. Recently, Tenneti et al. (2018) analyzed galaxies from the MassiveBlackII 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 masstolight 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 nonnegligible. 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 scaleinvariant at accelerations a < a_{0}, where a_{0} corresponds to g_{†} in Eq. (3) (Milgrom 2009). Thus, Eq. (3) is related to the interpolation function ν(g_{bar}∕a_{0}) 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_{0})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 nonzero 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. 
Maximum posterior parameters and reduced χ^{2} of individual rotation curve fits to the RAR.
References
 Battaglia, G., Fraternali, F., Oosterloo, T., & Sancisi, R. 2006, A&A, 447, 49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Begeman, K. G., Broeils, A. H., Sanders, R. H. 1991, MNRAS, 249, 523 [NASA ADS] [CrossRef] [Google Scholar]
 Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, S., & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton University Press) [Google Scholar]
 Blanchet, L., & Le Tiec, A. 2008, Phys. Rev. D, 78, 024031 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchet, L., & Le Tiec, A. 2009, Phys. Rev. D, 80, 023524 [NASA ADS] [CrossRef] [Google Scholar]
 Bosma, A. 1978, PhD Thesis, 1978, Groningen University, The Netherlands [Google Scholar]
 Bottema, R., & Pesta na J. L. G. 2015, MNRAS, 448, 2566 [NASA ADS] [CrossRef] [Google Scholar]
 Brada, R., & Milgrom, M. 1995, MNRAS, 276, 453 [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
 Burrage, C., Copeland, E. J., & Millington, P. 2017, Phys. Rev. D, 95, 064050 [NASA ADS] [CrossRef] [Google Scholar]
 Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709 [NASA ADS] [CrossRef] [Google Scholar]
 Chashchina, O., Foot, R., & Silagadze, Z. 2017, Phys. Rev. D, 95, 023009 [NASA ADS] [CrossRef] [Google Scholar]
 de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648 [NASA ADS] [CrossRef] [Google Scholar]
 Desmond, H. 2017a, MNRAS, 472, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Desmond, H. 2017b, MNRAS, 464, 4160 [CrossRef] [Google Scholar]
 Di Cintio, A., & Lelli, F. 2016, MNRAS, 456, L127 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Gentile, G., Famaey, B. & de Blok, W. J.G. 2011, A&A, 554, A125 [Google Scholar]
 Gentile, G., Józsa, G. I. G., Serra, P., et al. 2013, A&A, 527, A76 [Google Scholar]
 Katz, H., Lelli, F., McGaugh, S. S., et al. 2017, MNRAS, 466, 1648 [NASA ADS] [CrossRef] [Google Scholar]
 Keller, B. W., & Wadsley, J. W. 2017, ApJ, 835, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Khoury, J. 2015, Phys. Rev. D, 91, 024022 [NASA ADS] [CrossRef] [Google Scholar]
 Lelli, F., Verheijen, M., & Fraternali, F. 2014, A&A, 566, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016a, AJ, 152, 157 [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016b, ApJ, 816, L14 [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017a, ApJ, 836, 152 [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2017b, MNRAS, 468, L68 [NASA ADS] [CrossRef] [Google Scholar]
 Ludlow, A. D., BenítezLlambay, A., Schaller, M., et al. 2017, Phys. Rev. Lett., 118, 161103 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2004, ApJ, 609, 652 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2005, ApJ, 632, 859 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2014, Galaxies, 2, 601 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., & Schombert, J. M. 2014, AJ, 148, 77 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., & Schombert, J. M. 2015, ApJ, 802, 18 [CrossRef] [Google Scholar]
 McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok W. J. G. 2000, ApJ, 533, L99 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
 Meidt, S. E., Schinnerer, E., van de Ven, G., et al. 2014, ApJ, 788, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983, ApJ, 270, 371 [Google Scholar]
 Milgrom, M. 1994, Ann. Phys., 229, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2009, ApJ, 698, 1630 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2012, Phys. Rev. Lett., 109, 251103 [NASA ADS] [CrossRef] [Google Scholar]
 Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319 [Google Scholar]
 Navarro, J. F., BenitezLlambay, A., Fattahi, A. et al. 2017, MNRAS, 471, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Norris, M. A., Van de Ven, G., Schinnerer, E., et al. 2016, ApJ, 832, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Oman, K., Navarro, J. F., Fattahi, A. et al. 2015, MNRAS, 452, 3650 [CrossRef] [Google Scholar]
 Rubin, V. C., Thonnard, N., & Ford, Jr. W. K. 1978, ApJ, 225, L107 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & McGaugh, S. S. 2002, ARA&A, 40, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Schombert, J.,& McGaugh, S. S. 2014, PASA, 31, 36 [NASA ADS] [Google Scholar]
 Swaters, R. A., Sancisi, R., van Albada, T. S., & van der Hulst, J. M. 2009, A&A, 493, 871 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Swaters, R. A., Sancisi, R., van der Hulst, J. M., & van Albada, T. S. 2012, MNRAS, 425, 2299 [NASA ADS] [CrossRef] [Google Scholar]
 Tenneti, A., Mao, Y.Y., Croft, R. A. C., et al. 2018, MNRAS, 474, 3125 [NASA ADS] [CrossRef] [Google Scholar]
 Tully, R. B., & Fisher, J. R. 1977, A&A, 54, 661 [NASA ADS] [Google Scholar]
 Verheijen, M. A. W. 2001, ApJ, 563, 694 [NASA ADS] [CrossRef] [Google Scholar]
 Verlinde, E. P. 2017, SciPost Phys., 2, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, H., & Li, B. 2010, ApJ, 712, 130 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Maximum posterior parameters and reduced χ^{2} of individual rotation curve fits to the RAR.
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 dashdotted 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 gasdominated 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 lowquality 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 colorscale. 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 TullyFisher 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(M_{b}) = 4log(V _{f}) + log[X∕(a_{0}G)]. 

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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.