Testing Dark Matter and Modifications to Gravity using Local Milky Way Observables

TL;DR

This paper tests modified gravity models against Milky Way observables and finds that models mimicking MOND struggle to explain both rotation curves and vertical stellar motions, favoring dark matter explanations.

Contribution

It introduces a framework to evaluate scalar amplification gravity models using local Milky Way data, challenging MOND-like theories.

Findings

Models with scalar amplification are disfavored by Milky Way data.

Milky Way observations require substantial radial acceleration amplification.

MOND-like models cannot simultaneously explain rotation and vertical stellar motions.

Abstract

Galactic rotation curves are often considered the first robust evidence for the existence of dark matter. However, even in the presence of a dark matter halo, other galactic-scale observations, such as the Baryonic Tully-Fisher Relation and the Radial Acceleration Relation, remain challenging to explain. This has motivated long-distance, infrared modifications to gravity as an alternative to the dark matter hypothesis as well as various DM theories with similar phenomenology. In general, the standard lore has been that any model that reduces to the phenomenology of MOdified Newtonian Dynamics (MOND) on galactic scales explains essentially all galaxy-scale observables. We present a framework to test precisely this statement using local Milky Way observables, including the vertical acceleration field, the rotation curve, the baryonic surface density, and the stellar disk profile. We focus on models that predict scalar amplifications of gravity, i.e., models that increase the magnitude but do not change the direction of the gravitational acceleration. We find that models of this type are disfavored relative to a simple dark matter halo model because the Milky Way data requires a substantial amplification of the radial acceleration with little amplification of the vertical acceleration. We conclude that models which result in a MOND-like force struggle to simultaneously explain both the rotational velocity and vertical motion of nearby stars in the Milky Way.

Figures

Tables

Table 1. Table 1: Details of interpolation functions used in this study. Naming conventions, functional forms, priors and the values of Δ Δ \Delta BIC are given. The prior for the RAR function is taken directly from the results of fitting the function to the RAR in LABEL:McGaugh:2016leg , while the priors for the Simple and Standard functions are taken to have the same central value for a 0 subscript 𝑎 0 a_{0} , but with a somewhat larger variance to account for the fact that this value of a 0 subscript 𝑎 0 a_{0} is obtained by fitting with a different function. Note that here the notation LOGNORMAL ​ ( m , v ) LOGNORMAL 𝑚 𝑣 \mathrm{LOGNORMAL}(m,v) denotes a lognormal distribution with mean m 𝑚 m and variance v 𝑣 v . Values of Δ Δ \Delta BIC are given for a bulge mass M ∗ , bulge , obs = 1.50 ± 0.38 × 10 10 ​ M ⊙ subscript 𝑀 bulge obs plus-or-minus 1.50 0.38 superscript 10 10 subscript 𝑀 direct-product M_{*,{\rm bulge,obs}}=1.50\pm 0.38\times 10^{10}M_{\odot} .

Naming Convention	Functional Form	Prior for Scan	$\Delta$BIC
Taylor Expansion	$\nu \left(a_N\right)={\nu }_0+{\nu }_1{a}_N$	$\nu \left(a_N\right)>1\text{ or }1.3$	4.1 or 7.5
RAR McGaugh et al. (2016)	$\nu \left(a_N\right)={\left(1-e^{-\sqrt{a_N/a_0}}\right)}^{-1}$	$a_0=\mathrm{LOGNORMAL}(1.20,{0.24}^2)$	10.4
Simple Hees et al. (2016); Famaey and McGaugh (2012)	$\nu \left(a_N\right)=\frac{1}{2}\left(1+\sqrt{1+\frac{4}{a_N/a_0}}\right)$	$a_0=\mathrm{LOGNORMAL}(1.2,{0.4}^2)$	9.6
Standard Hees et al. (2016); Famaey and McGaugh (2012)	$\nu \left(a_N\right)=\sqrt{\frac{1}{2}\left(1+\sqrt{1+{\left(\frac{2}{a_N/a_0}\right)}^2}\right)}$	$a_0=\mathrm{LOGNORMAL}(1.2,{0.4}^2)$	4.8

Table 2. Table 2: Free parameters used in the benchmark Bayesian likelihood analysis along with the associated prior range for each and the 16-50-84 th th {}^{\text{th}} percentiles of the marginalized posterior. From left to right: normalization and scale height of i th superscript 𝑖 th i^{\text{th}} tracer population, normalization and scale length of stellar disk, bulge mass, normalization of gas disk, dark matter density normalization and inner slope, and interpolation function parameters for the ML models. The notation for the prior functions is: [ a , b ] 𝑎 𝑏 [a,b] denotes a flat prior between a 𝑎 a and b 𝑏 b ; LN ​ ( m , v ) LN 𝑚 𝑣 \mathrm{LN}(m,v) denotes a lognormal distribution with mean m 𝑚 m and variance v 𝑣 v .

	${\tilde{n}}_i$	$h_i$	${\tilde{\rho }}_{\ast }$	$h_{\ast ,R}$	$M_{\ast ,\mathrm{bulge}}$	${\tilde{\rho }}_g$	${\tilde{\rho }}_{\mathrm{DM}}$	$\alpha$	${\nu }_0$	$-{\nu }_1$	$a_0$
Unit	${10}^{-3}{\mathrm{pc}}^{-3}$	kpc	$M_{\odot }{\mathrm{pc}}^{-3}$	kpc	${10}^{10}{M}_{\odot }$	$M_{\odot }{\mathrm{pc}}^{-3}$	$M_{\odot }{\mathrm{pc}}^{-3}$	–	–	${10}^{10}{s}^2{m}^{-1}$	${10}^{-10}m{s}^{-2}$
Prior	[${10}^{-2}$, 1]	[${10}^{-3}$, 8]	[0, ${10}^2$]	[${10}^{-3}$, 8]	[0, ${10}^2$]	[0, ${10}^2$]	[0, ${10}^2$]	[0, ${10}^2$]	[0, ${10}^2$]	[0, ${10}^2$]	LN$(1.20,{0.24}^4)$
Posterior (ML)	${34.7}_{-4.8}^{+5.7}$ ${4.43}_{-0.41}^{+0.48}$ ${2.07}_{-0.12}^{+0.13}$	${0.26}_{-0.01}^{+0.01}$ ${0.46}_{-0.03}^{+0.03}$ ${0.87}_{-0.05}^{+0.06}$	${1.37}_{-0.95}^{+1.87}$	${2.41}_{-0.48}^{+1.26}$	${4.29}_{-2.37}^{+2.36}$	${0.25}_{-0.11}^{+0.14}$	–	–	${1.44}_{-0.08}^{+0.13}$	${0.21}_{-0.03}^{+0.05}$	–
Posterior (DM)	${33.6}_{-4.6}^{+5.4}$ ${4.35}_{-0.39}^{+0.44}$ ${2.06}_{-0.12}^{+0.13}$	${0.26}_{-0.01}^{+0.01}$ ${0.46}_{-0.03}^{+0.03}$ ${0.87}_{-0.05}^{+0.06}$	${1.21}_{-0.81}^{+1.17}$	${2.54}_{-0.45}^{+1.29}$	${2.87}_{-1.82}^{+2.12}$	${0.23}_{-0.10}^{+0.10}$	${2.16}_{-1.06}^{+1.12}$	${0.40}_{-0.30}^{+0.54}$	–	–	–
Posterior (ML RAR)	${39.8}_{-5.4}^{+6.5}$ ${4.8}_{-0.43}^{+0.48}$ ${2.14}_{-0.12}^{+0.13}$	${0.25}_{-0.01}^{+0.01}$ ${0.44}_{-0.02}^{+0.02}$ ${0.82}_{-0.05}^{+0.05}$	${0.56}_{-0.31}^{+0.52}$	${3.12}_{-0.62}^{+1.33}$	${4.10}_{-1.46}^{+1.05}$	${0.16}_{-0.05}^{+0.07}$	–	–	–	–	${0.91}_{-0.13}^{+0.14}$

Equations

$$a=\{\begin{array}{ll}{a}_{\mathrm{N}}& a\gg a_0\\ \sqrt{a_0{a}_{\mathrm{N}}}& a\ll a_0,\end{array}$$

$$\mathbf{a}=\nu \left(\frac{a_{\mathrm{N}}}{a_0}\right){\mathbf{a}}_{\mathrm{N}}.$$

$$\nu \left(\frac{a_{\mathrm{N}}}{a_0}\right)$$

$$\mathcal{L}({\mathbf{\theta }}_{\mathcal{M}})\propto \exp \left[-\frac{1}{2}\sum _{j=1}^N{\left(\frac{X_{j,\mathrm{obs}}-X_j({\mathbf{\theta }}_{\mathcal{M}})}{\delta X_{j,\mathrm{obs}}}\right)}^2\right]$$

$${\rho }_{\mathrm{B}}={\rho }_{\ast ,\mathrm{bulge}}+{\rho }_{\ast ,\mathrm{disk}}+{\rho }_{g,\mathrm{disk}}.$$

$${\rho }_{j,\mathrm{disk}}(R,z)={\tilde{\rho }}_j\exp \left(-R/h_{j,R}-|z|/h_{j,z}\right),$$

$${\rho }_{\ast ,\mathrm{bulge}}(r)=\frac{M_{\ast ,\mathrm{bulge}}}{2\pi }\frac{r_{\ast ,\mathrm{bulge}}}{r}\frac{1}{(r+r_{\ast ,\mathrm{bulge}}{)}^3},$$

$$\nabla \cdot \mathbf{a}=-4\pi G({\rho }_{\mathrm{B}}+{\rho }_{\mathrm{DM}}),$$

$${\rho }_{\mathrm{DM}}(r)=\frac{{\tilde{\rho }}_{\mathrm{DM}}}{{\left(r/r_s\right)}^{\alpha }{\left(1+r/r_s\right)}^{3-\alpha }},$$

$$\mathbf{a}=\left({\nu }_0+{\nu }_1{a}_{\mathrm{N}}\right){\mathbf{a}}_{\mathrm{N}},$$

$$\nabla \cdot {\mathbf{a}}_{\mathrm{N}}=-4\pi G{\rho }_{\mathrm{B}}.$$

$${\Sigma }_j^{z_{\max }}=2{\int }_0^{z_{\max }}{\rho }_j\left(R_0,z^{\prime }\right)d{z}^{\prime }.$$

$$v_c(R)={\sqrt{R\cdot a(R)}|}_{z=0}.$$

$$n_i(z)={\tilde{n}}_i\text{exp}(-\left|z\right|/h_i),$$

$${\sigma }_{i,z}(z{)}^2$$

$${\mathbf{\theta }}_{\mathrm{ML}}$$

$${\mathbf{\theta }}_{\mathrm{DM}}$$

$${\text{X}}_{\mathrm{obs}}=$$

$$M_{\ast ,\mathrm{disk}}=\frac{2\pi h_{\ast ,R}^2{\Sigma }_{\ast ,\mathrm{obs}}^{z_{\max }}\exp (R_0/h_{\ast ,R})}{1-\exp (-z_{\max }/h_{\ast ,z})},$$

$$h_{\ast ,R,\mathrm{obs}}=2.6\pm 0.5\text{leavevmode}\text{kpc}\text{and}{h}_{\ast ,z,\mathrm{obs}}=310\pm 50\text{leavevmode}\text{pc}.$$

$$\Delta \text{BIC}={\text{BIC}}_{\mathrm{ML}}-{\text{BIC}}_{\mathrm{DM}}.$$

$$r^2\to R^2+{\left(\frac{z}{c/a}\right)}^2.$$

$$\mathbf{a}=-\nabla \Phi ,{\nabla }^2\Phi =-\nabla \cdot \left(\nu \left(\frac{a_N}{a_0}\right){\mathbf{a}}_{\mathbf{N}}\right).$$

$$\mathbf{a}=\nu \left(\frac{a_{\mathrm{N}}}{a_0}\right){\mathbf{a}}_{\mathrm{N}}+\mathbf{S}.$$

$$\mathbf{S}=0\leftrightarrow \nabla |\nabla {\Phi }_{\mathrm{N}}|\times \nabla {\Phi }_{\mathrm{N}}=0.$$

$$a_{\mathrm{loc}}=\frac{v_c^2}{R_0}\approx 2\cdot 10^{-10}{\text{ m/s}}^2.$$

$${\Phi }_{\mathrm{N}}=-a_{\mathrm{N},\mathrm{BG}}z-\frac{G{m}_{\mathrm{pert}}}{r},$$

$$\Phi =-a_{\mathrm{BG}}z-\frac{G{m}_{\mathrm{pert}}}{r}\nu \left(\frac{a_{\mathrm{N},\mathrm{BG}}}{a_0}\right)\left[1+\frac{{\gamma }_{\mathrm{BG}}}{2}{\sin }^2\theta \right],$$

$$D\gg \sqrt{\nu \left(\frac{a_{\mathrm{N},\mathrm{BG}}}{a_0}\right)\frac{G{m}_{\mathrm{pert}}}{a_{\mathrm{loc}}}}.$$

$$D\gg 0.1\text{ kpc}\cdot {\left[\nu \left(\frac{a_{\mathrm{N},\mathrm{BG}}}{a_0}\right)\cdot \frac{m_{\mathrm{pert}}}{10^7{M}_{\odot }}\cdot \frac{2\cdot 10^{-10}{\text{ m/s}}^2}{a_{\mathrm{loc}}}\right]}^{1/2},$$

Full text

††preprint: PUPT XXXX

I Introduction

In the late twentieth century, the discovery that rotation curves flatten at large radii Rubin et al. (1980); Bosma (1981) revolutionized our understanding of galactic dynamics. This singular observation is often considered the first definitive evidence for the missing mass problem. The standard resolution assumes that a (mostly) collision-less dark matter (DM) halo provides the required increase in acceleration. However, the lack of direct evidence for DM to date, together with numerous galactic-scale observations, provides an opening for other resolutions of this puzzle, such as long-distance, infrared modifications to gravity as well as DM models that resemble such modifications on galactic scales. In this paper, we present a general framework to test the standard DM solution against such alternatives. One classic example of such an alternative to DM is MOdified Newtonian Dynamics (MOND) Milgrom (1983a, b, c). As a first application of our approach, we demonstrate the tension of models that reduce to a MOND-like force with local Galactic observables, thereby challenging them precisely at the scale where they are designed to be successful.

Many aspects of the well-known observation of flat rotation curves have yet to be understood—a statement that remains true even when considering DM. Examples include the Baryonic Tully-Fisher relation Tully and Fisher (1977), the Radial Acceleration Relation (RAR) McGaugh et al. (2016, 2000); Verheijen (2001) and additional correlations between baryons and observed accelerations on a galaxy-by-galaxy basis Sancisi (2003); McGaugh (2014); Lelli et al. (2017). While such correlations are somewhat challenging to explain within a DM paradigm, some theories of modified gravity provide a natural explanation because they avoid the introduction of additional matter. While these models often suffer from lack of a consistent short-distance, ultraviolet counterpart Bruneton and Esposito-Farese (2007) and struggle to explain cosmological observations Milgrom (2015), they are usually considered successful on galactic scales. MOND is probably the most well-known of these alternative explanations; it posits a modification to Newtonian gravity that becomes sizable only at extremely small accelerations. Even if not a fundamental theory of nature, MONDian phenomenology is extremely successful on galactic scales. It is therefore imaginable that the true theory, be it DM, a modification to gravity or some hybrid, resembles MOND on these small scales. One example which has recently attracted attention is Superfluid DM Berezhiani and Khoury (2015); Famaey et al. (2018); Berezhiani et al. (2018). This model results in an emergent MOND-like force on galactic scales while remaining in the DM phase on larger scales. Other novel theories which resemble MOND include TeVeS Bekenstein (2004), Emergent Gravity Verlinde (2011, 2017) and MOG Moffat (2006). Throughout this paper we will refer to any theory which reduces to MOND-like phenomenology on galactic scales as a “ML” model, even when the theory is not strictly a modification of gravity. While this study is mostly agnostic to the details of any ML model, the specific case of Superfluid DM will be the focus of an upcoming paper.

In this work, we use local stellar dynamics to test DM against a general class of ML models. The Milky Way (MW) rotation curve is essentially constant at the Solar radius Eilers et al. (2018), where the naive, baryon-only prediction falls off quite rapidly. As a result, local stellar measurements should be sensitive to the theory that governs the dynamics of the outer regions of the Galaxy Famaey and Binney (2005); Stubbs and Garg (2005); Nipoti et al. (2007); Bienayme et al. (2009); Famaey and McGaugh (2012); Bienaymé et al. (2014); Margalit and Shaviv (2016). Although DM and ML models can both successfully explain the observed flatness of the MW rotation curve, they do so through different mechanisms that enhance acceleration. In a spherical DM halo, the increase in radial acceleration is due to additional mass surrounding the Galactic disk and is independent of the baryonic distribution. By contrast, in the ML framework considered in this study, the additional acceleration is determined purely by the baryonic distribution and is often in the form of an approximate scalar enhancement to the Newtonian acceleration, i.e., an enhancement that changes the magnitude, but not the direction, of the Newtonian acceleration. In particular, for a spherically symmetric DM distribution, the net vector enhancement of the acceleration always points towards the Galactic Center, while in a ML scenario, the vector enhancement is in the direction of the Newtonian acceleration that is expected from baryons only.

The ability of our analysis to differentiate between these models comes from the fact that we constrain the radial and vertical (perpendicular to the midplane) accelerations simultaneously. Importantly, measurements of the baryonic profile of the MW and of local velocity dispersions Zhang et al. (2013) indicate that very little enhancement is required of the acceleration in the vertical direction. On the other hand, measurements of the MW rotation curve Eilers et al. (2018) require a much larger enhancement in the radial direction. Thus, a theory that treats both directions in an equivalent manner cannot easily explain all observations simultaneously, unless an anomalously large amount of baryonic mass is present at the center of the Galaxy. We find that local observations point towards a model that enhances the radial acceleration without significantly affecting the vertical acceleration. This provides a powerful handle to distinguish DM and any model that predicts a scalar-enhanced acceleration.

For the current study, we test ML forces using the local profile of the rotation curve, the local baryonic surface density, the vertical acceleration field within $\sim 1$ kpc of the Galactic midplane, and the shape of the MW stellar and gas profiles. These are standard observables that are often used to constrain the MW potential and to estimate the local DM density—see Refs. Bovy and Tremaine (2012); Garbari et al. (2012); Smith et al. (2012); Zhang et al. (2013); Bienaymé et al. (2014); Silverwood et al. (2016); Xia et al. (2016); Sivertsson et al. (2018) for the most recent studies, and LABEL:Read:2014qva for a review of the literature. We have performed a Bayesian likelihood analysis using these local MW observables to compare ML models and DM. In both cases, we model the MW baryonic density profile as a stellar disk, a gas disk, and a stellar bulge. We then fit these models to observations while marginalizing over uncertainties in the parameters. Tests of MOND using Galactic dynamics exist in the literature Iocco et al. (2015); Famaey and Binney (2005); Stubbs and Garg (2005); Nipoti et al. (2007); Bienayme et al. (2009); Famaey and McGaugh (2012); Bienaymé et al. (2014); Loebman et al. (2014); Margalit and Shaviv (2016), some of which have noted its effect on local vertical dynamics. However, a fully self-consistent study marginalizing over uncertainties in the baryonic distribution and directly comparing the results to a DM model has not been performed. Our analysis is also novel in that it is independent of any specific formulation of the ML model. Furthermore, our formalism has the potential to be extended to other scenarios.

We find that a theory with a ML force prefers a baryonic density profile that is in tension with known measurements of the stellar disk scale radius and the stellar bulge mass Bland-Hawthorn and Gerhard (2016); Licquia and Newman (2015). In particular, ML models require an anomalously small disk scale radius and/or an anomalously large bulge mass to reproduce local Galactic observables. On the other hand, DM is able to reproduce MW observables with parameters that are more consistent with the literature. Interestingly, the goodness of fit is slightly improved for a marginally prolate halo, i.e., one that introduces even less vertical acceleration than a spherical halo and is thus even further from the ML prediction.

Our results can be interpreted as model-independent in the sense that they need not rely on the details of the high and low-acceleration regimes in the Galaxy, but solely on the dynamics inferred from local measurements. When remaining completely agnostic to these details, the preference that we find for DM over ML models is positive, but not strong.111This study identifies observations where improvements in uncertainties can further strengthen the conclusions. On the other hand we find a strong preference for DM over a ML force for specific forms of the function which interpolates the high and low-acceleration regimes. One example of an interpolating function that is strongly constrained by our analysis is motivated by the recent study of rotation curves of $\sim 150$ galaxies in the SPARC database Lelli et al. (2016). Using this galaxy sample, Ref. McGaugh et al. (2016) finds the RAR, a universal relation between the observed radial acceleration and that expected due to baryons alone. The small scatter of the RAR can be interpreted as a manifestation of a ML force at play. Our result enables us to test the consistency of such an interpretation with Galactic observations. We find that the interpolation function suggested by Ref. McGaugh et al. (2016) is inconsistent with the observed vertical acceleration of disk stars near the Solar position. This does not mean, however, that the RAR is in tension with our results more generally.

The paper is organized as follows. Sec. II describes the general framework proposed by this study, focusing on the differences between a DM model and a ML model. Sec. III describes the Bayesian likelihood analysis that we perform. The modeling of the baryonic components are described in detail here, as well as the observational constraints that are used in the study. Sec. IV presents the results of the analysis, explicitly demonstrating that the best-fit stellar parameters in the ML model are in tension with observations, and discussing the systematic uncertainties that affect these results in detail. We conclude in Sec. V. The Appendix supplements the discussion on non-linear effects in MOND and ML models and also includes further details on the likelihood analysis.

II Dark Matter vs. MOND-Like Models

DM and ML models exhibit extremely different phenomenologies, even on galactic scales, due to their effects on the radial and vertical accelerations of tracer stars. In this section, we quantify this statement in further detail and point out specifically where the tension for ML models appears. Furthermore, we present details of our framework, which allows for strong distinguishing power between the different cases in a model-independent fashion.

As discussed above, the constraining power arises from the fact that for a well-motivated baryonic profile, little additional vertical acceleration is required to explain observations, while a relatively large increase in the radial acceleration must be invoked. The vertical acceleration is inferred from observations of the velocity dispersions of stars, while the radial acceleration is inferred from measurements of the local circular velocity of the MW.

DM is able to accommodate this requirement since an approximately spherical halo induces an acceleration pointing towards its center. Importantly, the net enhancement of the acceleration over the Newtonian contribution is just $\boldsymbol{a}_{\rm DM}$, the acceleration due to the DM component only, which is independent of the baryonic density profile. For a spherical DM profile, the local enhancement in cylindrical $(\hat{R},\hat{z})$ coordinates is $\boldsymbol{a}_{\rm DM}\approx-GM(R_{0})/R_{0}^{2}\times(1,z/R_{0})$, where $M(R_{0})$ is the enclosed DM mass, $z$ is the height above the midplane and $R_{0}$ is the radius of the Solar position. Thus, most of the local enhancement occurs in the radial direction with the vertical enhancement suppressed by $\sim z/R_{0}$.

The situation is quite different in a ML scenario, where the galactic dynamics are driven solely by the baryonic distribution. For many such models, the response to matter is highly non-linear, making a prediction of the dynamics at a given point within the galaxy non-trivial to calculate. However, as we demonstrate in this work, one need only characterize some general properties of the gravitational response to matter to provide discriminating power between ML models and DM. The particular characteristics that are important to classify are: (1) the parametric functional relationship between the local dynamical acceleration, $\boldsymbol{a}$, and the baryonic matter distribution, $\rho_{\rm B}$; (2) the tensor structure of this function; and (3) the degree to which this function varies in the region of interest within the MW.

For example, in many formulations of ML models, the observed acceleration, $\boldsymbol{a}$, depends only on $\rho_{\rm B}$ via the Newtonian acceleration, $\boldsymbol{a}_{\rm N}$. The asymptotic behavior is designed to reproduce the observed flatness of rotation curves in galaxies and is determined by

[TABLE]

where $a_{0}$ is a constant acceleration scale that sets the deviation from Newtonian gravity. In general, the solution for $\boldsymbol{a}$ is model-dependent and difficult to obtain. However, under certain conditions (see App. A), the dynamics reduces to the following form:

[TABLE]

The function $\nu(a_{\rm N}/a_{0})$ is known as the interpolation function and satisfies the asymptotic conditions of Eq. (1). Thus, for the case of a ML model, Eq. (2) is the parametric functional relationship between the dynamical acceleration and the baryonic distribution and it manifests as a scalar enhancement to the Newtonian acceleration. Importantly, in vanilla MOND, the specific form of this function is arbitrary, while in other models, e.g. Superfluid DM, it is set by the parameters of the model and/or the baryonic profile Berezhiani et al. (2018).

Finally, if one is only interested in fitting a ML model to Galactic observables in a localized region near the Sun, then the Newtonian acceleration, $a_{\rm N}$, does not vary too much from some particular reference value, $a_{\rm N,ref}$. Therefore, one can expand $\nu(a_{\rm N}/a_{0})$ in a Taylor series, giving

[TABLE]

to first order, where we have parametrized in terms of the constants $\nu_{0}\equiv\nu(a_{\rm N,ref}/a_{0})-\nu^{\prime}(a_{\rm N,ref}/a_{0})\cdot a_{\rm N,ref}/a_{0}$ and $\nu_{1}\equiv\nu^{\prime}(a_{\rm N,ref}/a_{0})/a_{0}$.222Note that in this parametrization, one can only constrain the value of $\nu_{1}$, i.e., $\nu^{\prime}(a_{\rm N,ref}/a_{0})/a_{0}$, and not the scale $a_{0}$ itself. In doing so, the arbitrary function in Eq. (2) is reduced to two constants, which can be treated as free parameters of the model. This expansion allows one to test any ML model and present results in a model-independent manner.

Given the phenomenology of each of these types of models, one can now ask which is better suited to simultaneously fit observations of the local radial and vertical accelerations. Measurements of these values are often inferred indirectly from the local circular velocity, $v_{c}(R_{0})$, which correlates with the radial acceleration through $a_{R}=v_{c}^{2}/R_{0}$, and from the vertical velocity dispersion, $\sigma_{z}(z)$, which correlates with the vertical acceleration through the Jeans equation Binney and Tremaine (2008). Fig. 1 presents the different predictions of these values for DM versus a ML model for $v_{c}(R_{0})$ and for $\sigma_{z}(z_{\rm ref})$ at a reference value of $z_{\rm ref}=1048$ pc.333The measurement of $\sigma_{z}(z_{\rm ref})$ used in this plot corresponds to the measurement of the metal-poor tracer population at $z_{\rm ref}=1048$ pc taken from Ref. Zhang et al. (2013), however the qualitative results shown in the figure are independent of these choices. For the purposes of illustration, we take a fixed baryonic profile that is consistent with measurements of the stellar disk and bulge as well as the gas disk. The black point indicates the measured values from Refs. Eilers et al. (2018) and Zhang et al. (2013).

The solid blue curve in the figure indicates the prediction for a DM model consisting of a spherical Navarro-Frenk-White (NFW) halo Navarro et al. (1996), varying over its normalization. Clearly, given this baryonic profile, the DM model easily accommodates the measured $v_{c}(R_{0})$ and $\sigma_{z}(z_{\rm ref})$ simultaneously. The dashed blue curve indicates the same correlation for a halo that is slightly prolate. We take the halo model of Ref. Bovy et al. (2016) with their best-fit axis ratio of $c/a=1.05$. Compared to the spherical case, this model enhances $v_{c}(R_{0})$ for a given $\sigma_{z}(z_{\rm ref})$ and is therefore slightly closer to the central values of the measurements.

The red solid curve in the figure is the correlation of these values for a ML model using Eq. (2) and the approximation for $\nu(a_{\rm N}/a_{0})$ in Eq. (3), varying over $\nu_{0}+\nu_{1}a_{\rm N}$. This curve shows that such models are unable to fit the measured values simultaneously. Specifically, in order to fit to $\sigma_{z}(z_{\rm ref})$, they under-predict $v_{c}(R_{0})$. One way around this is to invoke a different baryonic potential with much more mass at the center of the Galaxy. However, as is shown in the rest of this study, this is in conflict with measurements of the MW bulge mass and stellar disk scale radius.

The example presented in this section is meant for simple illustration of the primary tension underlying the ML and DM scenarios. To draw a more robust conclusion, one should marginalize over the baryonic uncertainties, which we describe in full in the next few sections.

Lastly, we comment that for ML models, the general relation between observed and Newtonian acceleration is slightly more complex than Eq. (2). Additionally, the non-linear structure of these theories can potentially cause further complications regarding the validity of the assumed smoothness of the baryonic profile. However, as is detailed in App. A, under certain circumstances, deviations from these assumptions are expected to be negligible. We find this to be the case for local measurements and thus for all results in this study. We stress that our analysis can be extended to more general frameworks; however, in this paper, we will restrict our focus to ML scenarios that satisfy Eq. (2).

III Methodology

This section describes the Bayesian likelihood analysis that we perform to test the consistency of the DM and ML models with local MW observables. For a given model $\mathcal{M}=\rm{DM},\rm{ML}$ with parameters $\boldsymbol{\theta}_{\mathcal{M}}$, we can predict the values of different MW observables $\mathbf{X}(\boldsymbol{\theta}_{\mathcal{M}})=(X_{1}(\boldsymbol{\theta}_{\mathcal{M}}),...,X_{N}(\boldsymbol{\theta}_{\mathcal{M}}))$. The predicted values are then compared to the measured values $\mathbf{X}_{\rm{obs}}=(X_{1,\rm{obs}},...,X_{N,\rm{obs}})$ with uncertainties $\delta\mathbf{X}_{\rm{obs}}=(\delta X_{1,\rm{obs}},...,\delta X_{N,\rm{obs}})$. We define the likelihood function to be

[TABLE]

and seek to recover the posterior distributions of the model parameters using Bayesian inference.444Note that in this equation we have assumed that the measurements are uncorrelated. Since most of the observations are taken from entirely different data sets, the correlations are expected to be small. An analysis which accounts for the full covariance matrix is beyond the scope of this study.

In Sec. III.1, we define the models, in Sec. III.2 we describe the observables $\mathbf{X}_{\rm obs}$ that are used as constraints, and in Sec. III.3 we discuss the analysis procedure in detail.

III.1 Model Definition

In this section, we detail the parametrization of the DM and ML models. In both cases, the predicted acceleration depends sensitively on the MW baryonic mass density, $\rho_{\rm B}$, which we take to consist of a stellar bulge and a disk of stellar ($*$) and gaseous ($g$) components:

[TABLE]

For both the stellar and gaseous disks, we use a double-exponential density profile given by

[TABLE]

where $h_{j,R}$ ($h_{j,z}$) is the scale length (height) of the disk, $\tilde{\rho}_{j}$ is the normalization, and the index $j=*,g$ corresponds to either the stellar or gas component. The variables $R$ and $z$ are Galactocentric cylindrical coordinates.

For computational simplicity, we do not separately model the thin and thick stellar disks, but rather treat these contributions together as a single exponential density profile with an effective scale height and radius. As the local density ratio of the thick to thin disk ranges from $\sim 1$–10% Gilmore and Reid (1983); Siegel et al. (2002); Juric et al. (2008), the properties of the effective disk are primarily set by the thin component. We note that the effective disk can also be considered to be the weighted contribution of individual mono-abundance populations Bovy et al. (2012); Bovy and Rix (2013), which yields roughly the same effective scale length at the Solar position compared to what we estimate from the thin/thick disk contributions.

In principle, there are 6 free parameters for the disk model; however, we restrict half of these values. For our benchmark analysis, we take $h_{g,z}=130$ pc and $h_{g,R}=2\,h_{*,R}$, following Refs. Binney and Tremaine (2008); Bovy and Rix (2013), and also fix the stellar scale height to $h_{*,z}=300$ pc, consistent with measurements from stellar counts Bland-Hawthorn and Gerhard (2016).

For the stellar bulge, we use a Hernquist density profile Hernquist (1990),

[TABLE]

where $r_{*,\rm bulge}$ is the bulge scale radius, and $M_{*,\rm bulge}$ is a normalization constant that sets the total mass of the bulge. Here, $r$ is the Galactocentric spherical radial coordinate. The scale radius is fixed to $r_{*,\rm bulge}=600$ pc, while the normalization is allowed to vary. While the MW bulge has more structure than captured by Eq. (7) (see LABEL:Bland-Hawthorn:2016aaa for a review), our analysis is not sensitive to these details as the observables of interest are measured beyond $R\gtrsim 5$ kpc.

Finally, we must specify the dynamics of the DM and ML models. For the DM scenario, the dynamical acceleration is given by the solution to Poisson’s equation,

[TABLE]

where $\rho_{\rm DM}$ is the DM mass density. In particular, we assume that the DM density profile follows a generalized NFW distribution Navarro et al. (1996),

[TABLE]

where $r_{s}$ is the scale radius, and $\alpha$ is the inner slope. The parameters $\tilde{\rho}_{\mathrm{DM}}$ and $\alpha$ are both free in our analysis procedure, while we fix $r_{s}=19$ kpc.

On the other hand, for the ML scenario in our benchmark analysis, the observed dynamical acceleration is given by

[TABLE]

where $\nu_{0}$ and $\nu_{1}$, defined in Eq. (3), parametrize the interpolation function and $\mathbf{a}_{\rm N}$ is the Newtonian acceleration, which satisfies

[TABLE]

We allow $\nu_{0}$ and $\nu_{1}$ to be free parameters in the analysis procedure; however, we restrict the model to ensure that the observed acceleration is always an enhancement of Newtonian gravity by forcing ${\nu_{0}+\nu_{1}a_{\rm N}>1}$ everywhere within the spatial region of interest.

In addition to this benchmark analysis, we also consider a number of specific interpolation functions, as well as the use of the interpolation function of Eq. (10) with the requirement that ${\nu_{0}+\nu_{1}a_{\rm N}>1.3}$. Details of all functions are presented in Table 1. For the interpolation functions which are not Taylor expansions, there is a single free parameter, $a_{0}$, for which we set priors as detailed in the table. The results of these analyses are far less general than the use of Eq. (10), but are more in line with specific formulations of MOND and ML models, including Superfluid DM. Generically, they force $\mathbf{a}>\mathbf{a}_{\rm N}$ at the Solar position, whereas with Eq. (10) and the requirement ${\nu_{0}+\nu_{1}a_{\rm N}>1}$, the model also allows $\mathbf{a}\approx\mathbf{a}_{\rm N}$ at the Solar position.

III.2 Constraints from Local Observations

Next, we discuss the values of the observables that are used in the likelihood analysis. The solar radius, $R_{0}$, feeds into the prediction for each of these observables.555In this work, we neglect the vertical offset of the Sun above the Galactic plane, i.e., we take $z_{0}=0$. This parameter is critical because any Galactic measurements based on angular size or velocity perceived from Earth depend on $R_{0}$. We take $R_{0}=8.122$ kpc as the fiducial value, consistent with the observation of the orbit of the star S2 around the massive black hole candidate Sgr A* Abuter et al. (2018).

The first constraint arises from the local stellar and gas surface densities, as determined from photometric observations. By definition, the surface density depends on the distribution of the stellar or gas component as

[TABLE]

We adopt measurements of ${\Sigma_{g,\text{obs}}^{1.1}=12.6\pm 1.6\ \mathrm{M}_{\odot}\ \mathrm{pc}^{-2}}$ and ${\Sigma_{*,\text{obs}}^{1.1}=31.2\pm 1.6\ \mathrm{M}_{\odot}\ \mathrm{pc}^{-2}}$ at ${z_{\rm max}=1.1\ \mathrm{kpc}}$ as fiducial values McKee et al. (2015). For self-consistency, we only use surface densities obtained from direct photometric observations (as opposed to dynamical studies), which are equally valid for the DM and ML scenarios.

The second constraint comes from the value and slope of the rotation curve at the Solar radius. The circular velocity and its derivative with respect to $R$ are obtained from the predicted acceleration using

[TABLE]

We consider only the local circular velocity and slope of the rotation curve at $R=R_{0}$, taking $v_{c,\text{obs}}=229\pm 12\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ sec^{-1}$ and $\left(dv_{c}/dR\right)_{\rm obs}=-1.7\pm 0.47\leavevmode\nobreak\ \rm km\leavevmode\nobreak\ sec^{-1}\leavevmode\nobreak\ kpc^{-1}$ as fiducial values Eilers et al. (2018). Note that we have added in quadrature the statistical and systematic uncertainties quoted in LABEL:2018arXiv181009466E and, in particular, that the $\sim 5\%$ systematic uncertainties in $v_{c}$ dominate.

The final constraint comes from the observed number density and vertical velocity dispersions of three mono-abundance stellar populations at $R=R_{0}$, provided in LABEL:Zhang:2012rsb. These results were obtained using 9000 K-dwarfs in the SEGUE sub-survey of the Sloan Digital Sky Survey (SDSS) York et al. (2000); Abazajian et al. (2009); Aihara et al. (2011). The sample was categorized based on its iron fraction, [Fe/H], and $\alpha$-abundance, [$\alpha$/Fe], and divided into three tracer populations consisting of metal-rich, metal-intermediate, and metal-poor stars. For each population, indexed by $i$, number densities, $n_{i,{\rm obs}}(z_{k})$, and vertical velocity dispersions, $\sigma_{z,i,{\rm obs}}(z_{k})$, were obtained for several values of $z_{k}$ between $300\leavevmode\nobreak\ \mathrm{pc}$ and $1200\ \mathrm{pc}$.

Following LABEL:Zhang:2012rsb, we model the number densities, $n_{i}(z)$, for the $i^{\text{th}}$ tracer population as

[TABLE]

where $\tilde{n}_{i}$ and $h_{i}$ are additional parameters varied in the fit. The tracer stars are assumed to be in steady state and to be well-described by the Jeans equations Binney and Tremaine (2008). In this case, the vertical velocity dispersion is

[TABLE]

where $a_{z}$ is the predicted vertical acceleration. Although the integral in Eq. (15) extends to infinity, the integrand falls off exponentially in $z$ due to Eq. (14), and we are not sensitive to dynamics above $z\gtrsim h_{i}\sim 1\ \mathrm{kpc}$. Additionally, Eq. (15) ignores the contribution from the tilt term, which depends on the velocity dispersion in the $R-z$ plane. The contribution of this term is negligible, as estimated in LABEL:Zhang:2012rsb for the specific K-dwarf sample studied there, and can thus be safely neglected.

III.3 MCMC Analysis

To summarize, the complete set of model parameters is formed by the union of the set of SEGUE tracer population parameters, $\{\tilde{n}_{i},\,h_{i}\,|\,i\leavevmode\nobreak\ =\leavevmode\nobreak\ 1,2,3\}$, the baryon potential parameters, $\{\tilde{\rho}_{*},\ h_{*,R},\ h_{*,z},\ M_{*,\rm bulge},\ r_{*,\rm bulge},\ \tilde{\rho}_{g},\ h_{g,R},\ h_{g,z}\}$, and $\{\nu_{0},\nu_{1}\}$ (or $a_{0}$) for the ML model or $\{\tilde{\rho}_{\text{DM}},\alpha,r_{s}\}$ for the DM case. We set $r_{s}$, $r_{*,\rm bulge}$, $h_{*,z}$, and $h_{g,z}$ to constant values, and fix $h_{g,R}=2h_{*,R}$. It would be preferable to marginalize over all parameters in the scans, however this is computationally prohibitive. The effects of varying these fixed values are summarized in detail in Sec. IV.4. The free parameters in Eq. (4) are therefore

[TABLE]

(for scans with specific interpolation functions $\{\nu_{0},\ \nu_{1}\}\to a_{0}$) and the local observables are

[TABLE]

For a given model, $\mathcal{M}$, the predicted values of these observables, $\textbf{X}(\boldsymbol{\theta}_{\mathcal{M}})$, are obtained from Eqs. (12)-(15).

We use the Affine Invariant Markov Chain Monte Carlo (MCMC) implementation emcee Foreman-Mackey et al. (2013) to recover the posterior distributions for $\boldsymbol{\theta}_{\rm ML}$ and $\boldsymbol{\theta}_{\rm DM}$, as appropriate, using the likelihood defined in Eq. (4). For our benchmark analysis, we set flat linear priors on each parameter as summarized in Table 2 and additionally require that each walker step satisfy (i) $\nu_{0}+\nu_{1}a_{\rm N}>1$ and (ii) that the baryonic rotation curve peak below $R<5$ kpc.666This is motivated since the MW is known to have a maximal disk Sackett (1997), with the rotation curve peak occurring in a region that is baryon-dominated. For the additional analyses, we set either flat or log-normal priors as detailed in Table 1. The details of the MCMC analysis, including posterior distributions, two-parameter correlations, and convergence, are presented in App. B.

IV Results

IV.1 The Baryonic Density Profile

This section presents the results of the Bayesian parameter inference described in Sec. III. The median values and uncertainties of all the free parameters of our baseline scans (as well as the RAR function scan for comparison) are presented in Table 2. We start by exploring the preferred baryonic mass distributions for both the DM and ML scenarios. The results are summarized in Fig. 2, which shows the correlations between the stellar disk scale radius $h_{*,R}$, stellar disk mass $M_{*,{\rm disk}}$, and bulge mass $M_{*,{\rm bulge}}$. Note that, as will be discussed below, in the final results presented at the end of this section, the scale radius and bulge mass are set as constraints in the analysis. However, the correlations shown in the figure best illustrate the result of this study: the DM model is more consistent with observations than the ML model. The 68% and 95% containment regions are indicated for both the ML (red) and DM (blue) models.

The left panel of Fig. 2 shows the correlation of $h_{*,R}$ with $M_{*,{\rm disk}}$, for a specific value of the stellar scale height, $h_{*,z}=300$ pc. The correlation between these two parameters is mostly driven by the constraint on the local stellar surface density. Using Eqs. (6) and (12), it can easily be shown that

[TABLE]

where $z_{\rm max}=1.1$ kpc for the constraints used in this study. The resulting curve is indicated by the black dashed line; both the ML and DM best-fit regions roughly follow this trend.

The central panel of Fig. 2 illustrates the correlation of the stellar bulge mass, $M_{*,{\rm bulge}}$, with the stellar disk mass, $M_{*,{\rm disk}}$. This correlation is driven primarily by the value of the local circular velocity, which is roughly determined by the total enclosed baryonic mass. Therefore, in the limit of spherical symmetry and neglecting the gas disk, $M_{*,\rm disk}+M_{*,\rm bulge}$ should be constant to recover the observed circular velocity; thus, $M_{*,\rm disk}$ and $M_{*,\rm bulge}$ are negatively correlated. This argument can be confirmed more quantitatively for the case of cylindrical symmetry as follows. For each value of $M_{*,\rm disk}$, we compute the unique baryon profile that gives the observed circular velocity and stellar and gas surface densities, fixing baryonic parameters as in the baseline analysis and holding the DM or ML parameters fixed. Specifically, for the ML case, we hold $\nu_{\odot}\equiv\nu_{0}+\nu_{1}\cdot a_{\rm N}(R_{0})=1.3$,777In the ML case, the circular velocity is only sensitive to $\nu_{0}$ and $\nu_{1}$ in the form of $\nu_{\odot}$. approximately consistent with the median value from the scan (which is $\nu_{\odot}=1.1$). This is in line with the intuition described in the previous sections: the ML scan picks out an acceleration enhancement that is close to unity. The results are plotted as the dashed lines in the central panel of Fig. 2, and we see that they agree appreciably well with the intuition outlined above.

The right panel of Fig. 2 shows the correlation of $h_{*,R}$ with $M_{*,{\rm bulge}}$—a convolution of the results in the first two panels. As $h_{*,R}$ increases, the disk contribution to $a_{{\rm N},R}$ decreases, causing $M_{*,\rm bulge}$ to increase in order to maintain the constant $v_{c,\rm obs}$. The DM and ML regions both show this general trend. The dashed lines representing the DM and ML estimates are obtained following the same procedure we used to get the estimates in the central panel.

The open/filled black points in the right panel of Fig. 2 represent two different measurements of the bulge mass: $M_{*,{\rm bulge,obs}}=1.50\pm 0.38$ $M_{\odot}$ from microlensing observations Calchi Novati et al. (2008) (filled black circle) and $M_{*,{\rm bulge,obs}}=0.62\pm 0.31$ $M_{\odot}$ (open black circle) from photometric observations Lopez-Corredoira et al. (2007). Note that we have rescaled the latter to the value of $R_{0}$ used in this work Licquia and Newman (2015). In both cases, the disk scale length is taken to be $h_{*,R,{\rm obs}}=2.6\pm 0.5$ kpc. It is important to emphasize that we model the stellar disk using a single exponential profile, and thus $h_{*,R,\rm obs}$ is an effective value derived from a two-disk (thin and thick) formalism. We perform a Monte Carlo procedure that properly combines the weighted contributions from the thin and thick disks. For the thin disk, we adopt $h_{*,R,\rm obs}^{t}=2.6\pm 0.5$ kpc and $h_{*,z,\rm obs}^{t}=300\pm 50$ pc, and for the thick disk we adopt $h_{*,R,\rm obs}^{T}=2.0\pm 0.2$ kpc and $h_{*,z,\rm obs}^{T}=900\pm 180$ pc, as derived in LABEL:Bland-Hawthorn:2016aaa from a review of the literature. We obtain the following estimates for the effective disk parameters:

[TABLE]

We note that $h_{*,z,\rm obs}$ is consistent with, although not precisely equal to, the benchmark value of $h_{*,z}=300\ \text{pc}$.

As mentioned above, in our final results, we will fix the scale length and bulge mass directly in the likelihood analysis. However, for illustrative purposes and because of the spread of values found in the literature (see LABEL:Bland-Hawthorn:2016aaa and LABEL:Licquia:2014rsa for a summary), we have also considered the case where these parameters are unconstrained in order to understand the general preference of the DM and ML models. Thus, the results presented in Fig. 2 are quite general, and one can easily interpret them in the context of any specific measurement for the disk parameters. From this figure, we see that ML models prefer larger $M_{*,\rm bulge}$ and smaller $h_{*,R}$, as compared to DM, and have less overlap with the favored region indicated by the open/filled black points. In order to understand this behavior, it is enlightening to study posterior distributions of a number of ML parameters.

Fig. 3 shows the posterior distributions of $\nu(a_{\rm N}(R_{0}))$ for two ML models. The red region corresponds to the baseline analysis with full freedom of the interpolation function, while the grey region corresponds to the analysis which enforces the RAR interpolation function (see Table 1). The inset shows the prior and posterior distributions of $a_{0}$ for the RAR function. Each model attempts to reduce the acceleration enhancement as much as is allowed within the adopted prior range. For example, the red region is pushed to the minimum allowed value of $\nu(a_{\rm N}(R_{0}))$ as allowed by the priors on $\nu_{0}$ and $\nu_{1}$ (i.e., $\nu_{0}+\nu_{1}a_{\rm N}>1$). The gray region is also forced to small values of $\nu(a_{\rm N}(R_{0}))$, with the location of the peak set by the fact that the $a_{0}$ posterior is peaked towards its lowest allowed value. Because of the functional form of the interpolation function, as $a_{0}$ decreases, the local acceleration enhancement (the value of $\nu(a_{\rm N}(R_{0}))$) also decreases, i.e., accelerations become more Newtonian. The gray region is still peaked above the red because it is associated with a specific interpolation function and preferred value of $a_{0}$. This forces a model with the RAR interpolation function to significantly amplify local accelerations, making it challenging for this model to simultaneously explain radial and vertical accelerations. On the other hand, the Tayler expansion function is able to explain the vertical accelerations by choosing $\nu(a_{\rm N}(R_{0}))\approx 1$, but then requires anomalous amounts of baryonic mass at the center of the Galaxy to explain the radial accelerations.

Fig. 4 plots the discrepancy of each step in the MCMC sample chain against the measured values $h_{*,R,\rm obs}$ and $M_{*,{\rm bulge,obs}}$;888We compute this discrepancy assuming that both measurements are uncorrelated and normally distributed. The discrepancy is then $\sqrt{\left(\frac{h_{*,R}-h_{*,R,\rm obs}}{\delta h_{*,R,\rm obs}}\right)^{2}+\left(\frac{M_{*,\rm bulge}-M_{*,\rm bulge,obs}}{\delta M_{*,\rm bulge,obs}}\right)^{2}}$, where $\delta h_{*,R,\rm obs}\leavevmode\nobreak\ (\delta M_{*,\rm bulge,obs})$ is the measurement error of the scale length (bulge mass). the solid (dashed) lines correspond to the larger (smaller) bulge mass measurement. This result is shown only for the baseline DM and ML scans. The majority of the walker steps fall within $\sim 1\sigma$ of the $h_{*,R}$ observation for the DM model, while the corresponding distribution for the ML case is peaked around $\sim 2\sigma$ discrepancy. This figure is an alternative way to represent the results from the right panel of Fig. 2, and highlights the fact that the DM model is more consistent with observations of the baryonic profile.

We can also use the Bayesian Information Criterion (BIC) Kass and Raftery (1995) to compare the DM and ML models, defining

[TABLE]

The larger the value of $\Delta\text{BIC}$, the more strongly disfavored the ML model is relative to the DM model. To obtain these values, we run the scans imposing additional constraints on the parameters. If we take the scale length to be that of Eq. (19) and the bulge mass to be $1.50\pm 0.38\times 10^{10}$ ($0.62\pm 0.31\times 10^{10}$) $M_{\odot}$, then $\Delta\text{BIC}=4.1$ (4.5) for the baseline scan where the local interpolation function is constrained only to be larger than 1. These values suggest a positive, but not strong, preference for the DM model. Notice that the $\Delta\text{BIC}$ value increases when the bulge mass is smaller. We have also calculated the $\Delta\text{BIC}$ for scans with a number of interpolation functions as detailed in Table 1. The Taylor expansion function with $\nu(a_{\rm N})>1$ corresponds to the results given above. The other scans result in larger $\Delta\text{BIC}$ values (given for a bulge mass of $1.50\pm 0.38\times 10^{10}$ $M_{\odot}$) as follows. For the Taylor expansion function with $\nu(a_{\rm N})>1.3$, $\Delta\text{BIC}=7.5$. For the RAR function, $\Delta\text{BIC}=10.4$. For the Simple function, $\Delta\text{BIC}=9.6$. For the Standard function, $\Delta\text{BIC}=4.8$. All these results more strongly favor a DM model. The underlying reason being that these functions enforce a local acceleration enhancement substantially larger than unity and the tension described in the previous sections becomes more apparent.

Finally, in line with the reasoning described in Sec. II, we have performed a scan with a DM halo that is slightly prolate. In this case, the ratio of $a_{{\rm DM},z}/a_{{\rm DM},R}$ is even more suppressed relative to a ML model. We take the halo to have the functional form described in Ref. Bovy et al. (2016), which is equivalent to Eq. (9) with the replacement

[TABLE]

We take the best-fit axis ratio from that study, namely $c/a=1.05$. In this scan, we impose the constraints $M_{*,{\rm bulge}}=1.50\pm 0.38\times 10^{10}$ $M_{\odot}$ and $h_{*,R}=2.6\pm 0.5$ kpc. This model is preferred over a ML model with $\Delta\text{BIC}=4.2$. We have verified that allowing $c/a$ to scan to larger values increases this preference.

IV.2 Kinematic Comparison of the Models

As discussed in Sec. II, the primary tension between ML models and MW observables lies in the inability to simultaneously reproduce the radial and vertical accelerations of stars. As a test of this, we consider the correlation between the ratios $a_{z}/a_{\text{N},z}$ and $a_{R}/a_{\text{N},R}$, evaluated at $R=R_{0}$ and at $z=300$ and 1200 pc. Fig. 5 plots these two ratios for the models under consideration. The ML model directly relates each of these ratios to the interpolating function, $\nu(a_{\rm N}/a_{0})$, so they are always equivalent. In other words, in a ML model, the ratio of the observed and the Newtonian accelerations is independent of direction. The ML model prediction is simply the solid red line in Fig. 5.

The DM case is comparatively more interesting. The difference between the DM confidence regions and the ML expectation highlights the main kinematic difference between the two models, as anticipated in the discussion of Sec. II. If the DM best-fit region were to partially overlap with the ML line, then the models would be kinematically indistinguishable in a certain region of parameter space. We instead see that there is always a sizable difference between the two cases, for all values of $z$ between 300–1200 pc. We also note that since the slope of the DM correlation should be set by $z/R_{0}$, the DM and ML kinematics become more and more similar as $z$ is increased. This is because the disk potential looks more spherical as $z$ increases. In the limit of a spherically symmetric baryonic distribution, we expect ML models and DM to be indistinguishable, at least in the sense described in Sec. II.

IV.3 Dark Matter and MOND-like Parameters

In addition to allowing us to compare the DM and ML models, our Bayesian analysis framework allows us to recover the best-fit parameters specific to each. For example, the right panel of Fig. 6 shows the 68% and 95% containment region for the local DM density, $\rho_{\text{DM}}(R_{0})$, and the inner slope of the NFW profile, $\alpha$. Although $\rho_{\rm DM}(R_{0})$ is not an independent free parameter of our framework, it is easily computable from Eq. (9). We find that the marginalized value for the local DM density is $\rho_{\rm DM}(R_{0})=0.29\pm 0.06\ \mathrm{GeV}\ \mathrm{cm}^{-2}$, consistent with previous measurements summarized in LABEL:Read:2014qva. In particular, we can compare our results to those of LABEL:Zhang:2012rsb, whose K-dwarf number densities and vertical dispersions are used here (see Sec. III.2). The two studies differ in several important analysis details—for example, we vary over the full baryonic profile and assume an NFW distribution. However, the final results are similar, with LABEL:Zhang:2012rsb finding a local DM density of $0.25\pm 0.09$ GeV cm*-3*. We also constrain the fully marginalized value of the inner slope of the NFW profile to be $\alpha<1.1$ at 90% confidence. This result, which is consistent with previous dynamical studies Bovy and Rix (2013), is not constraining enough to distinguish between a cored ($\alpha\approx 0$) and cusped profile ($\alpha\approx 1$).

The left panel of Fig. 6 shows the analogous containment for the ML parameters, $\nu_{0}$ and $\nu_{1}$, in the red shaded region, for the baseline scan with $\nu(a_{\rm N})>1$. We find that the fully marginalized values are $\nu_{0}=1.4\pm 0.1$ and $\nu_{1}=-0.21^{+0.03}_{-0.05}\times 10^{10}\ \mathrm{s}^{2}\ \mathrm{m}^{-1}$. Because $\nu_{0}$ and $\nu_{1}$ are a model-independent parametrization of the interpolation function, these values can be related to any specific interpolation function. As an example, we consider the RAR interpolation function of LABEL:McGaugh:2016leg, which is given in Table 1 with $a_{0}=1.20\pm 0.24\times 10^{-10}\ \mathrm{m}\ \mathrm{s}^{-2}$ (note that we have now chosen the specific central value and errors of $a_{0}$ which is the best fit to the RAR data). In order to map the interpolation function onto the $\nu_{0}$ and $\nu_{1}$ parameters, the expansion parameter of Eq. (3), $a_{\rm N,ref}$, must be specified. Although in principle $a_{\rm N,ref}$ can be evaluated at any point within the local region of interest, it will differ very little from $a_{\rm N}(R_{0},z_{0})$; we therefore approximate $a_{\rm N,ref}\approx a_{\rm N}(R_{0},z_{0})$.999We have verified that the gray region in the left panel of Fig. 6 does not change appreciably when $a_{\rm N,ref}$ is evaluated within $|R-R_{0}|\lesssim 500\ \mathrm{pc}$ and $|z|\lesssim 2\ \mathrm{kpc}$. We then perform the following Monte Carlo procedure to estimate the values of $\nu_{0}$ and $\nu_{1}$ that correspond to the RAR interpolation function with this value of $a_{0}$: for each step in the MCMC sample chain for the ML model, we compute $a_{\rm N}(R_{0},z_{0})$, generate a random value of $a_{0}$ that is normally distributed with mean and standard deviation as specified by LABEL:McGaugh:2016leg, and then compute $\nu_{0}$ and $\nu_{1}$ using Eq. (3). The gray shaded regions in the left panel of Fig. 6 indicate 68% and 95% containment for the recovered values of $\nu_{0}$ and $\nu_{1}$. We note that the region of maximum posterior probability for $\nu_{0}$ and $\nu_{1}$ in our analysis is beyond the gray shaded countours, indicating that our analysis disfavors the interpolation function of LABEL:McGaugh:2016leg together with its best-fit values for $a_{0}$. We note, however, that this result is specific to this specific interpolation function. A different function could simultaneously fit the RAR and be consistent with our results.

We stress that our treatment of the interpolation function from LABEL:McGaugh:2016leg is purely demonstrative. Given any interpolation function and value of $a_{0}$, one can check its consistency against our results. To this end, we have included the fully marginalized posterior distribution for $a_{\rm N}(R_{0},z_{0})$ in the inset of the left panel of Fig. 6. The distribution is well-approximated by a Gaussian with mean $1.65\times 10^{-10}\ \mathrm{m}\ \mathrm{s}^{-2}$ and standard deviation $0.22\times 10^{-10}\ \mathrm{m}\ \mathrm{s}^{-2}$. Although we used the full MCMC sample chain to construct the gray region in Fig. 6, we have verified that approximating the $a_{\rm N}(R_{0},z_{0})$ distribution with the aforementioned Gaussian does not significantly affect the results and should serve as a sufficient estimator for one who seeks to perform this procedure without access to the complete MCMC sample chain.

IV.4 Systematic Checks

Lastly, we explore the effects of varying the baseline analysis by repeating the MCMC scan for different values of the fixed parameters. To assess the overall impact on the final conclusions, we impose the constraints $h_{*,R,{\rm obs}}=2.6\pm 0.5$ kpc and $M_{*,{\rm bulge,obs}}=1.50\pm 0.38\times 10^{10}$ $M_{\odot}$. We quote the $\Delta\text{BIC}$ as defined in Eq. (20), and compare to the baseline value of 4.1, for the following variations:

• •

fix the bulge scale radius to $r_{*,\rm bulge}=2$ kpc, rather than $600$ pc. ($\Delta\text{BIC}=3.7$)

• •

fix the NFW scale radius to $r_{s}=25$ kpc, rather than 19 kpc. ($\Delta\text{BIC}=3.9$)

• •

fix the gas scale length to $h_{g,R}=4.5$ kpc Kalberla (2003); Barros et al. (2016), rather than requiring $h_{g,R}=2\,h_{*,R}$ kpc ($\Delta\text{BIC}=4.6$). For comparison, the median values of the stellar scale length for the ML (DM) baseline study is $h_{*,R}=2.41\leavevmode\nobreak\ (2.54)$ kpc, corresponding to a median gas scale length of $h_{g,R}=4.82\leavevmode\nobreak\ (5.08)$ kpc.

• •

fix the gas scale height to $h_{g,z}=180$ pc, rather than 130 pc. ($\Delta\text{BIC}=4.2$)

• •

change the stellar disk scale height to $h_{*,z}=200,400,600$ pc, compared to 300 pc. ($\Delta\text{BIC}=4.1,3.5,2.9$ respectively)

• •

allow $h_{*,z}$ to vary as a free parameter of the scan and apply the constraint that $h_{*,z,{\rm obs}}=310\pm 50$ pc. ($\Delta\text{BIC}=4.1$)

• •

use logarithmic, rather than linear, priors for all parameters in the scan. ($\Delta\text{BIC}=3.4$)

In general, we find that the conclusions are not strongly sensitive to the variations listed above, with the results always providing positive evidence in favor of DM. Of the parameters that we varied over, the analysis is perhaps the most sensitive to the stellar disk scale height. To understand this behavior, we observe that increasing the disk scale height, while fixing the local stellar surface density, decreases the stellar density close to the midplane. As a result, the vertical velocity dispersion decreases as well, forcing a larger value of $\nu_{\odot}$ in order for a ML model to match observations. Because the radial acceleration obeys the scaling $a_{R}\simeq\nu_{\odot}a_{R,N}$, a larger $\nu_{\odot}$ is compensated by a smaller $a_{R,N}$ to reproduce the observed circular velocity. This, in turn, implies a smaller stellar bulge mass and/or larger stellar scale radius, in better agreement with observations.

As an additional variation, we replace the effective stellar disk with two double exponential profiles—one each for the thin and thick disk—to verify that the use of a more complicated density profile does not dramatically affect the results. For simplicity, we fix the thick disk parameters and the thin disk scale height using the best-fit values from the SDSS study in LABEL:Juric:2005zr, and we keep the scale length of the thin disk as a free parameter in the MCMC scan. The gas disk parameters are treated in the same way as in the baseline study. In this case, we find $\Delta\text{BIC}=4.1$.

An important assumption entering into the Jeans analysis described in Sec. III.2—specifically Eq. (15)—is that the tracer stars are in steady state. This assumption is challenged by recent evidence of waves in nearby disk stars, possibly caused by the passage of a dwarf galaxy such as Sagittarius through the disk Widrow et al. (2012); Williams et al. (2013); Gomez et al. (2013); Antoja et al. (2018); Bennett and Bovy (2019). Such perturbations can lead to $\mathcal{O}(10\%)$ changes to the vertical force Banik et al. (2017). We have performed systematic checks to see how robust our results are to the uncertainties in the constraints used. For example, we repeated the baseline MCMC scan, (arbitrarily) doubling the errors on the tracer number density and vertical dispersion data from LABEL:Zhang:2012rsb. While we still find positive evidence in favor of DM ($\Delta\text{BIC}=3.1$), the preference is smaller. This is expected because larger errors would reduce the tension between the ML prediction and the vertical acceleration data. As an additional test, we (artificially) double the errors on the surface density constraint taken from LABEL:McKee:2015hwa; in this case, $\Delta\text{BIC}=4.1$, and is unchanged.

V Conclusions

In this work, we tested the consistency of DM and models which resemble MOND on galactic scales with local MW observables. Our goal was to compare the models on galactic scales, where ML models have had their greatest success. As a concrete example, we focused on scalar enhancements to Newtonian gravity, where the absolute magnitude, but not the direction, of the Newtonian acceleration is altered. In this case, the direction of the total acceleration vector is set by the baryonic mass distribution. This is fundamentally different than the DM expectation; in the presence of a spherical halo, the total acceleration vector necessarily points towards its center.

Such phenomenology is generically predicted by models which reduce to a ML force on galactic scales. One example is Superfluid DM Berezhiani and Khoury (2015); Famaey et al. (2018); Berezhiani et al. (2018), where the ML force arises as an emergent force in a superfluid phase of DM. In that model, there are additional complications arising from the existence of a small DM density together with the ML force. While we postpone the full treatment of these models to future work, we expect a qualitative similarity to the current study where we consider the case of a ML force alone. In general, we find that scalar enhancements to gravity over-predict the vertical acceleration of nearby stars when trying to simultaneously explain the observed radial acceleration. This provides a crucial handle for testing these models on galactic scales.

We have performed a Bayesian likelihood analysis to compare DM to ML models. The results depend sensitively on the baryonic mass distribution, which we describe using a stellar disk, a gas disk, and a stellar bulge. We varied over the parameters of the baryonic profile, imposing minimally restrictive priors. For the DM scenario, the halo distribution was assumed to follow a spherical NFW profile. For the ML scenario, we focused on the class of models where the baryonic acceleration is enhanced by a scalar function. For our baseline analysis, we parametrized this scalar interpolation function in a model-independent manner by considering its linear approximation near the Solar position. This approach allowed us to evaluate the preference for DM, independent of the details of the interpolation function. The results are generic and conservative, as the parametrization has total freedom to enhance accelerations as much or as little as required by the data (including for example no enhancement at all). In this case, we find positive, but not strong, preference for DM relative to ML models. However, when we compare DM to ML models with specific interpolation functions, we find strong tension with data.

The free model parameters were constrained using local measurements of the baryonic surface density, rotational velocity, and vertical dispersion of tracer stars. We note that, on general grounds, the non-linear nature of ML models might result in an acceleration field that is not parallel to the Newtonian prediction. However, as is demonstrated in App. A, we expect any deviations from alignment with the Newtonian acceleration to be small.

Our findings can be summarized as follows. In comparison to DM, ML models typically prefer a smaller scale length for the stellar disk and a larger stellar bulge mass, in tension with current measurements. The underlying reason for this is that MW data seem to prefer a model which is able to simultaneously enhance radial acceleration with little enhancement to vertical accelerations. ML models are generically unable to provide such phenomenology unless they invoke anomalously large amounts of baryonic mass at the Galactic Center. For our baseline analysis and taking an observed scale length of $h_{*,R,{\rm obs}}=2.6\pm 0.5$ kpc and bulge mass of $M_{*,{\rm bulge},{\rm obs}}=1.50\pm 0.38\times 10^{10}$ ($0.62\pm 0.31\times 10^{10}$) $M_{\odot}$, we find that $\Delta\text{BIC}=4.1$ (4.5), favoring the DM model. These conclusions are robust to systematic variations of the analysis. We also find that a slightly prolate DM halo is also consistent with the data. When the analysis is done with any specific interpolation function, forcing acceleration enhancements to be considerably larger than unity, the $\Delta\text{BIC}$ grows to values in the range 4.8–10.4, depending on the specific functional form and its prior.

The tension that we have observed for ML models versus DM will be further clarified as measurements of the baryonic density profile continue to improve. As we have seen, the disk scale length and bulge mass play a particularly important role in defining this tension, so reducing the measurement uncertainty on these parameters has the potential to strengthen the conclusions. Additionally, we note that the measurements of the MW rotation curve are systematics-dominated and that improvements in the measurement of the local circular velocity would also strengthen the power of the analysis by tightening the contour bands in Fig. 2. As a simple demonstration of this, we have rerun the analysis using only statistical errors on these quantities, ignoring the quoted systematic uncertainties. This increases the preference for the DM model, with $\Delta\rm{BIC}=7.1$.

This initial study does not take full advantage of the Gaia dataset Gaia Collaboration et al. (2018). Indeed, the only input that we have used from Gaia is the recent update to the rotation curve from LABEL:2018arXiv181009466E. As the characterization of local MW dynamics continues to improve with Gaia, one can revisit the analysis proposed in this work to obtain an even more definitive test of the models.

VI Acknowledgements

We thank J. Bovy, J. Herzog-Arbeitman, D. Hogg, G. Kribs, S. McGaugh, H. Walter-Rix, H. Verlinde, and L. Zhang for useful conversations. This research made use of the Astropy Astropy Collaboration et al. (2013), IPython Pérez and Granger (2007), matplotlib Hunter (2007), numpy Oliphant (2015), galpy Bovy (2015), corner Foreman-Mackey (2016), and emcee Foreman-Mackey et al. (2013) software packages. ML is supported by the DOE under contract DESC0007968 and the Cottrell Scholar Program through the Research Corporation for Science Advancement. MM is supported by the DOE under contract DESC0007968. NJO is supported by the Azrieli Foundation Fellows program.

Appendix A Non-Linear Effects in MOND and MOND-like models

A.1 The Divergenceless Field

Any analysis of the type presented in this study critically assumes some correlation between the Newtonian acceleration that is expected from the baryonic distribution and the observed acceleration. Because ML models typically invoke a non-linear response to the distribution of baryonic matter, any simple relation between Newtonian and observed acceleration may hold only approximately. If this is the case, then one must first ensure that any corrections do not significantly affect the analysis. For our analysis, we have assumed Eq. (2) to be true, and we now describe the possible corrections to this and why we believe them to be small.

Eq. (2) is in fact only a special case of the correlation between the observed and Newtonian acceleration, and is known to violate standard conservation laws Milgrom (1983b). To circumvent this, it is therefore required to define a scalar field whose gradient is the acceleration field, in analogy with Newtonian gravity. A simple realization is provided in Ref. Milgrom (2010) as

[TABLE]

When the Newtonian acceleration is a function of a single coordinate, it can be shown that the above equation reduces to Eq. (2) using the divergence theorem Brada and Milgrom (1995); Bekenstein and Milgrom (1984). More generally, however, Eq. (22) holds up to a divergence-less vector field, $\boldsymbol{S}\equiv\boldsymbol{\nabla}\times\boldsymbol{h}$, such that101010Strictly speaking, Eqs. (22,23) correspond to a version of MOND known as Quasi-Linear MOND (QuMOND) Milgrom (2010). However, both the standard formulation of MOND and QuMOND are equally motivated, and they reduce to each other exactly in the limit where $\boldsymbol{S}=0$.

[TABLE]

The MW baryonic profile is disk-like, and even in the limit of a perfect disk, two coordinates are required to describe it (namely $R$ and $z$). Thus, in principle, $\boldsymbol{S}$ need not be zero and other considerations must be used to ensure that any effects arising from a non-zero $\boldsymbol{S}$ field are negligibly small.

We are interested in ensuring that we can safely assume that $\boldsymbol{S}\ll\nu\left(a_{\rm N}/a_{0}\right)\boldsymbol{a}_{\rm N}$ in the region of study used in this work. Following the arguments of Ref. Brada and Milgrom (1995), since $\boldsymbol{S}$ is both divergenceless and irrotational, the following self-consistency condition relates the requirement of vanishing $\boldsymbol{S}$ with a requirement on the Newtonian potential:

[TABLE]

This equation is satisfied when $\boldsymbol{\nabla}\Phi_{\rm N}=0$ or when $\boldsymbol{\nabla}|\boldsymbol{\nabla}\Phi_{\rm N}|$ is parallel to $\boldsymbol{\nabla}\Phi_{\rm N}$, and can be verified explicitly for any Newtonian potential. Specifically, the same study found this requirement to hold extremely well for any baryonic profile with a radial exponentially decreasing surface density profile. We have repeated the same analysis assuming a baryonic profile that consists of a stellar disk ($M_{*,{\rm disk}}=3.5\pm 1\times 10^{10}$ $M_{\odot}$, $h_{*,R}=2.6\pm 0.5$ kpc, $h_{*,z}=300\pm 50$ pc), a gas disk ($M_{g,{\rm disk}}=1\pm 0.4\times 10^{10}$ $M_{\odot}$, $h_{g,R}=5.2\pm 1$ kpc, $h_{g,z}=130\pm 50$ pc) and a Hernquist stellar bulge ($M_{*,{\rm bulge}}=1\pm 0.4\times 10^{10}$ $M_{\odot}$), finding that any deviation from Eq. (24) cannot be distinguished to within measurement errors of $|\boldsymbol{\nabla}\Phi_{\rm N}|$. The values of the baryonic profile parameters used here are motivated by measurements of the bulge mass Licquia and Newman (2015), the stellar Bland-Hawthorn and Gerhard (2016) and gas Binney and Tremaine (2008); Bovy and Rix (2013) disk parameters ($M_{g,{\rm disk}}$ is the total mass of the gas disk) as well as local surface density measurements McKee et al. (2015). We have also varied over these parameters, and find that this result is not sensitive to any $\mathcal{O}(1)$ changes to their values. The advantage of this type of consideration is that it does not require any knowledge of $\nu(a_{\rm N}/a_{0})$ and is therefore robust. The disadvantage is that it is not straightforward to quantify the size of $\boldsymbol{S}$ given that Eq. (24) is only approximately satisfied.

A.2 Mass Perturbations

An additional complication arises from the fact that the MW baryonic density profile is not perfectly smooth. Because of the non-linear nature of MOND and ML models, small perturbations to a smooth matter distribution can potentially cause extremely large effects. This is related to the previously studied External Field Effect Bekenstein and Milgrom (1984); Milgrom (1983b); Famaey and McGaugh (2012); Haghi et al. (2009); Hees et al. (2016). Thus, if one is to model the MW with a smooth cylindrical density profile, it is of the utmost importance to ensure that small perturbations do not cause large effects in the region of study.

Near the Sun’s position in the MW disk, there are a number of known stellar overdensities Juric et al. (2008). These regions explicitly break assumptions regarding cylindrical symmetry and the smoothness of the baryonic density distribution. In Newtonian gravity, the acceleration field is dictated by the Poisson equation, whose linearity enables one to very simply determine when a perturbation can be neglected. The non-linear nature of ML models requires a more subtle treatment since, as will be shown below, the effects of a perturbation can be felt over larger distances.

To quantify this effect, the above situation can be modeled as follows. Consider a point mass, $m_{\rm pert}$, at a position within the disk, but sufficiently far from the Galactic Center such that the background Galactic acceleration field (in the vicinity of the point mass), $\boldsymbol{a}_{\rm BG}$, can be treated as constant. The point mass creates its own acceleration field which depends on $\boldsymbol{a}_{\rm BG}$ because of the non-linear behavior of the model. We are interested in quantifying the size of this field and comparing it to the local observed acceleration, $\boldsymbol{a}_{\rm loc}$, which can be directly inferred observationally, e.g., from the local value of the rotation curve,

[TABLE]

Setting the coordinate origin to the position of the point mass and aligning the $\hat{z}$ direction with $-\boldsymbol{a}_{\rm BG}$ (see Fig. 7), the Newtonian potential around the point mass can be written as

[TABLE]

where $a_{\rm N,BG}$ is the Newtonian background field that relates to the observed background field through $\nu(a_{\rm N,BG}/a_{0})\,a_{\rm N,BG}=a_{\rm BG}$. To leading order in $m_{\rm pert}$, the perturbed field (background and perturbation) is given by

[TABLE]

where $\gamma_{\rm BG}\equiv\partial\log\nu/\partial\log{[a_{\rm N}/a_{0}]}$ evaluated at $a_{\rm N,BG}/a_{0}$ Milgrom (2010). The result of Eq. (27) allows one to determine the maximal range of influence of the perturbation by demanding that its effect be much smaller than $a_{\rm loc}$. This requires that the distance, $D$, from the perturbation satisfies

[TABLE]

Note that $D$ is always larger than the equivalent distance for the Newtonian case where $\nu=1$.

There are a number of known stellar overdensities within $\mathcal{O}(1)$ kpc of the Sun. The most prominent of these are localized at $(R,Z)\sim(6.5,1.5)\,\rm{kpc}$ and $\sim(9.5,0.8)\,\rm kpc$, with densities of approximately twice that of their average surroundings and sizes of order kpc3—see, e.g., Fig. 27 of LABEL:Juric:2005zr. Comparing to the smooth density profile model found by these authors at the locations of these perturbations, we estimate the total mass of each overdensity to be $m_{\rm pert}\lesssim 10^{7}\,M_{\odot}$. To safely neglect the effects of these overdensities, the distance from the mass perturbation must be larger than

[TABLE]

which is automatically satisfied for any interpolation function whose value in the vicinity of the perturbation is $\nu\approx 1$.

Appendix B MCMC details and convergence

We run emcee with 200 walkers for 200,000 steps per walker, discarding the first 40,000 steps as burn-in. Based on visual inspection of the parameter trace plots, we expect this choice of burn-in to be highly conservative. We note that the limiting factor in burn-in time is the parameter $h_{*,R}$, which gradually diverges and eventually saturates its prior range. A corner plot of all the free parameters and their posterior distributions and two-parameter correlations can be seen for the ML model in Fig. 8 and for the DM model in Fig. 9. We assess the convergence of the post burn-in chain by computing the autocorrelation time for each parameter Foreman-Mackey et al. (2013); Goodman and Weare (2010). We estimate that the autocorrelation length for each parameter is $\lesssim$ 10,000 steps, which suggests that our sample chain is long enough to be a sufficient proxy for the true posterior distribution. In Fig. 10, we show the estimated autocorrelation length, $\hat{\tau}$, as a function of the number of post burnin-in steps, $N$, as a visual indicator of the convergence for both the ML and DM baseline scans. We repeat similar analyses for the additional emcee scans that we run with additional measurements $h_{*,R,\rm obs}$ and $M_{*,\rm bulge,obs}$ in the likelihood; in general, we find that these scans converge much faster.