Comparing dark matter and MOND hyphotheses from the distribution function of A, F, early-G stars in the solar neighbourhood

Abstract

Dark matter is hypothetical matter believed to address the missing mass problem in galaxies. However, alternative theories, such as Modified Newtonian Dynamics (MOND), have been notably successful in explaining the missing mass problem in various astrophysical systems. The vertical distribution function of stars in the solar neighbourhood serves as a proxy to constrain galactic dynamics in accordance to its contents. We employ both the vertical positional and velocity distribution of stars in cylindrical coordinates with a radius of 150 pc and a half-height of 200 pc from the galactic plane. Our tracers consist of main-sequence A, F, and early-G stars from the GAIA, RAVE, APOGEE, GALAH, and LAMOST catalogues. We attempt to solve the missing mass in the solar neighbourhood, interpreting it as either dark matter or MOND. Subsequently, we compare both hypotheses newtonian gravity with dark matter and MOND, using the Bayes factor (BF) to determine which one is more favoured by the data. We found that the inferred dark matter in the solar neighbourhood is in range of $\sim(0.01\textup{--}0.07)\,\textrm{M}_{\odot}\,\textrm{pc}^{-3}$. We also determine that the MOND hypothesis’s acceleration parameter $a_{0}$ is $(1.26\pm 0.13)\times 10^{-10}\,\textrm{m}\,\textrm{s}^{-2}$ for simple interpolating function. The average of bayes factor for all tracers between the two hypotheses is $\log\textrm{BF}\sim 0.1$, meaning no strong evidence in favour of either the dark matter or MOND hypotheses.

Keywords: Dark Matter – Gravitation – Galaxy: kinematics and dynamics

1   INTRODUCTION

The study of dark matter constitutes a crucial topic in modern astrophysics. The discovery of flat rotation curves at large radii in spiral galaxies offered a glimpse into the existence of dark matter (Rubin et al. 1980). The local density of dark matter, typically observed a few hundred parsecs from the Sun, holds significant importance in constraining the nature of dark matter within the Milky Way. This value then imparts information about the shape of the Milky Way’s dark matter halo in the midplane. Furthermore, it serves as a constraint in galaxy formation models and cosmology (Read 2014), and even extends its importance to alternative theories of gravity (Milgrom 2001; Nipoti et al. 2007).

One of the first documented explorations of the missing mass problem was conducted by Oort (1932), who examined the vertical force on stars in the vicinity of the Sun. He observed that the vertical force exceeded the expected value based on visible matter alone. A similar problem was also identified by Zwicky (1933) through the observation of the mass-to-light ratio in the Coma cluster. Zwicky found that the mass-to-light ratio derived from dynamical means significantly exceeded that derived from the luminous matter alone.

The mainstream explanation for the missing mass problem is the existence of non-baryonic dark matter. This hypothetical matter interacts solely through gravity with baryonic matter, resulting in an additional gravitational force among the stars. The adoption of dark matter particles in various astrophysical systems has proven highly successful in numerous areas. Examples include explaining the rotation curve of spiral galaxies (Rubin et al. 1980), ensuring the stability of galactic discs (Ostriker & Peebles 1973), addressing mass deficiencies in ultrafaint dwarf galaxies (Lin & Ishak 2016), and contributing to our understanding of cosmological evolutions (Davis 2000).

The local dark matter density is typically estimated by measuring the vertical force on stars near the Sun. This is commonly accomplished by measuring the vertical velocity, a practice that had been studied since Oort (1932). Read (2014) provides an extensive review of the methods used to measure local dark matter. The estimated value of the local dark matter density usually hovers around $\sim 0.01$ M${}_{\odot}$ pc${}^{-3}$, although some studies show higher values, reaching up to $\sim 0.1$ M${}_{\odot}$ pc${}^{-3}$ (Bahcall 1984a). Fig. 1 displays the values of local dark matter density from various studies.

Figure 1: The estimated local dark matter density from a few studies in the last decade, including our own. The points are from Moni Bidin et al. (2012; MB12), Garbari et al. (2011; G12), Zhang et al. (2013; Z13), Bienaymè et al. (2014; B14), Xia et al. (2016; X16), Sivertsson et al. (2018; S18), and Buch et al. (2019; B19) for type A, F, and early G stars.

Alternative studies suggest that the apparent effect of missing mass on the vertical (and radial) force may not be attributed to dark matter itself but rather to modifications in the gravitational force law. This is typically achieved by adjusting the Newtonian gravitational force law, introducing slight changes to the effect of gravity at large distances. One classic example of such a modification is the Modified Newtonian Dyanimcs (MOND) theory, proposed by Milgrom (1983). Subsequently, extensive studies have been conducted to replace dark matter with modified gravity theories, yielding varying degrees of success (Famaey & McGaugh 2012).

The MOND theory has proven highly successful in elucidating the dynamics of stars in spiral galaxies (Milgrom 1994). This theory mimics the effects attributed to dark matter, but notably, it does not necessitate the presence of additional matter. Nevertheless, there are some differing predictions between the MOND theory and the dark matter theory. Lisanti et al. (2019) attempt to exploit these differences, particularly in the ratio of the vertical force to the radial force in the Milky Way. Their findings reveal that the MOND theory predicts a distinct ratio compared to the dark matter theory, providing a means to test which theory aligns more closely with the observed data.

Another method employed to test the hypotheses of dark matter and MOND involves comparing the centripetal acceleration of wide binary stars with their acceleration calculated using classical Newtonian gravity. Hernandez (2023) discovered anomalous behavior in the low acceleration regime, favoring MOND. Chae (2023) even identified evidence at a $10\sigma$ level in favor of MOND. Conversely, Pittordis & Sutherland (2023) found that standard gravity is more preferred. A more rigorous Bayesian statistical analysis conducted by Banik et al. (2023) revealed results consistent with Newtonian dynamics, ruling out MOND with a confidence level of $16\sigma$. This complex array of findings indicates the need for further studies to rigorously test these two competing hypotheses.

In this paper, we will explore the influence of matter distribution (with or whitout dark matter) on generating gravitational potential and shaping the vertical distribution of stars near the Sun. We assume that the stars are in dynamical equilibrium and employ the collisionless Boltzmann equation to model the distribution of stars (Holmberg & Flynn 2004). The vertical density profile of tracers and their vertical velocity dispersion are measured. Subsequently, we attempt to fit the vertical distribution using mass models that assume isothermal components in dynamical equilibrium. This analysis is conducted for both dark matter and MOND hypotheses. The comparison of both best-fittings is then performed by evaluating their Bayes factor.

The paper is organized as follows: In Section 2, we provide a detailed description of the data collection and the subsequent data cleaning process. Section 3 outlines the mass model, theoretical background, and the two hypotheses. Section 4 explains the methodology employed for data fitting using Markov Chain Monte Carlo (MCMC). It includes an explanation of the likelihood function, the prior, and the numerical method used for solving the gravitational potential and calculating the Bayes factor. Section 5 presents a the analysis of the results, discussing the dark matter density from the dark matter hypothesis, the acceleration parameter from the MOND hypothesis, and a comparison of the two hypotheses using Bayes factors. Finally, Section 6 summarizes and concludes the findings of our study. Additionally, Section A contains the derivation of the likelihood function for inferring number density.

2   DATA

The data utilised in this study consist of the three-dimensional positions and kinematics of stars around the Sun. These data are compiled from Gaia DR3, RAVE 6, APOGEE-2 (SDSS), GALAH DR3, and LAMOST DR7 catalogues. They are collected either by downloading the bulk data directly from their repositories or queried using astroquery through virtual observatory services.

The Gaia catalogue serves as the primary source, providing three-dimensional positions ($\alpha,\delta,\varpi$) and kinematics ($\mu_{\alpha},\mu_{\delta},v_{r}$) of stars. Other catalogues are employed to fill in missing radial velocity data from Gaia DR3. Table 1 displays the distribution of radial velocity data from each catalogue. It is evident that the contribution of radial velocity data from catalogues other than Gaia is not significant, rendering the subsequent analysis unaffected by data collected solely from Gaia DR3.

Table 1: Combined contribution of radial velocity data from Gaia, RAVE, APOGEE, GALAH, and LAMOST catalogues.

Catalogue	#
GALAH	78,055
APOGEE	370,523
RAVE	27,826
LAMOST	442,140
Gaia	23,191,519

Additional information from 2MASS photometry is also utilised to obtain the colours $J$ and $K_{s}$ of the stars. The colour $J-K_{s}$ is used for selecting the main sequence cut, following a methodology similar to that employed in Bovy’s work (2017a). Infrared photometry is chosen because it is less susceptible to the effects of interstellar extinction. Hence, for this study, there is no need to correct for the extinction effect. The 2MASS All-sky Point Source Catalogue (PSC) encompasses photometry data for 470,992,970 stars in the near-infrared wavelength ($J$, $H$, and $K_{s}$ bands).

2.1 Sample Selection

Firstly, all stars from Gaia DR3 are selected across all directions (all sky). Subsequently, we crossmatch these stars with the 2MASS PSC to acquire the colours $J$ and $K_{s}$. This crossmatching process utilises the 2MASS-Gaia joined catalogue, a compilation previously conducted by the Gaia team (Marrese, P. M. et al. 2019) and accessible from the Gaia archive. Both catalogues are outer joined, combining entries whether there are matching stars from 2MASS to Gaia or not. The resulting catalogue, named Gaia+2MASS, encompasses astrometry, kinematics, and photometry data for 470,992,970 stars.

The remaining catalogues are obtained by downloading the data from their respective archives. Fortunately, for each catalogue, there exist a pre-existing crossmatch with Gaia DR3 conducted by the Gaia team. These crossmatches are established by simply matching the IDs from the Gaia+2MASS catalogue with the corresponding join table catalogue. Specifically, there are 450,978 crossmatched stars from RAVE, 573,733 stars from APOGEE, 588,061 stars from GALAH, and 2,730,472 stars from LAMOST. It is required for the crossmatched catalogue to include radial velocity data for the stars.

Further filtering is conducted to exclude stars with poor-quality photometry, astrometry, and outliers. Adopting the criteria outlined in Kim & Lépine’s work (2021), the following selection criteria are applied to the stars:

• Parallax $\varpi>0$: Only stars with a positive parallax are selected to avoid the complexities associated with distance inference for stars with negative parallax.

• Parallax error $\sigma_{\varpi}/\varpi<0.15$: Stars with good-quality parallax data, characterized by a relative error less than 15 per cent, are chosen.

• $-3.0<G_{\textup{BP}}-G_{\textup{RP}}<6.0$ and $3.0<G<21.0$: This criterion aim to eliminate outliers in the colour-magnitude diagram.

• $\sigma\left(F_{\textup{BP}}\right)/F_{\textup{BP}}<0.10$ and $\sigma\left(F_{\textup{RP}}\right)/F_{\textup{RP}}<0.10$: Only stars with good-quality photometry in $G_{\textup{BP}}$ and $G_{\textup{RP}}$ bands, characterized by relative errors less than 10 per cent, are selected.

• ruwe $<1.4$: Only stars with good astrometry, defined by a renormalized unit weight error (ruwe) less than 1.4, are chosen.

• $1.0+0.015\left(G_{\textup{BP}}-G_{\textup{RP}}\right)^{2}<E<1.3+0.06\left(G_{\textup{BP}}-G_{\textup{RP}}\right)^{2}$: This criterion is applied to select stars with a good colour excess $E$.

Additional filtering is carried out by incorporating the stars’ infrared photometry from the 2MASS PSC to obtain the highest-quality infrared data. Following Bovy’s methodology (2017a), the criteria for selecting the best infrared data are:

• $2<J<13.5$: Establishing a lower bound to eliminate outliers and an upper bound to ensure greater than 99 per cent completeness.

• ph_qual = A_A: Selecting high-quality photometry based on the photometric quality flag.

• rd_flg = 1 or 3: Verifying data quality using the read flag, specifically selecting stars with flag values of 1 or 3.

• use_src $=1$: Choosing stars that satisfy the criteria for good 2MASS data.

All the aforementioned criteria are applied to the Gaia+2MASS catalogue, yielding 29,144,488 stars with Gaia ID, referred to as the pre-combined catalogue. Additionally, 50,250,399 stars are identified with 2MASS infrared photometry only, serving as the reference dataset.

In the final step, the pre-combined catalogue is crossmatched with the remaining catalogues using the Gaia DR3 source ID. This crossmatched catalogue is denoted as the combined catalogue, encompassing numerous radial velocity data from various catalogues. The radial velocity data employed in this study corresponds to the radial velocity with the smallest error from the respective original catalogue. Stars lacking radial velocity data are excluded from the combined catalogue, resulting in a dataset comprising 24,110,063 stars with complete astrometry ($\alpha,\delta,\varpi$), kinematics ($\mu_{\alpha},\mu_{\delta},v_{r}$), and photometry ($J,K_{s}$) data.

2.2 Midplane Selection and Main Sequence Cut

Inferences regarding the stellar dynamics in the Milky Way due to matter distribution are drawn by analyzing the positions and motions of appropriate stellar populations, referred to as tracers. For this study, we specifically focus on stars in the solar neighborhood, defining them as stars within a cylindrical volume in galactice coordinate with a vertical distance of $|z|<200$ pc from the midplane and a radius of $150$ pc. Main sequence stars are chosen and categorized into multiple groups based on their colour $J-K_{s}$. Following the methodology used by Bovy (2017a), we employ the mean dwarf stellar locus from Pecaut & Mamajek (2013) and introduce a reasonable thickness to the absolute magnitude $M_{J}$. Utilising a similar approach, we define the main sequence by tracing the main branch of the Hertzsprung-Russell diagram (HRD) of nearby stars, establishing upper and lower boundaries in absolute magnitude as displayed by the shaded area in Fig. 2.

Figure 2: The Hertzsprung-Russell (HR) diagram depicts stars with stellar type A to early-G in the solar neighborhood, with the shaded area denoting the main sequence cut.

The stellar populations under consideration in this study must have achieved dynamical equilibrium within the Galactic disc in recent times. Following the recommendation of Buch et al. (2019), we opt for stars characterized by lower scale height and a shorter equilibration timescale, specifically focusing on younger A stars up to early G stars. Consequently, we choose the interval between A0 stars and G3 stars and divide the group into 15 subgroups, as detailed in Table 2. The color division for each group is carefully selected to ensure a minimum of 1000 stars in the midplane, guaranteeing robust statistical analyses for each subgroup. The spectral type of each group is determined from the color $J-K_{s}$ using the relation provided by Pecaut & Mamajek (2013).

Table 2: Tracer groups of stars utilised in the study, categorized by colour. The midplane column indicates the number of stars within $|z|<50$ pc.

$J-K_{s}$	Spectral Type	$|z|<200$ pc	Midplane
$[-0.019$, $0.124]$	A0 - A8	2422	1006
$[0.124$, $0.188]$	A9 - F1	2719	1014
$[0.188$, $0.219]$	F2 - F4	2867	1027
$[0.219$, $0.241]$	F5	2919	1014
$[0.241$, $0.259]$	F6	2970	1033
$[0.259$, $0.275]$	F7	3291	1049
$[0.275$, $0.288]$	F8	3165	1032
$[0.288$, $0.300]$	F9	3315	1038
$[0.300$, $0.312]$	F9.5	3523	1060
$[0.312$, $0.323]$	F9.5	3546	1001
$[0.323$, $0.333]$	G0	3493	1007
$[0.333$, $0.343]$	G1	3580	1028
$[0.343$, $0.353]$	G1	3803	1067
$[0.353$, $0.362]$	G2	3517	1044
$[0.362$, $0.376]$	G2 - G3	3460	1529

2.3 Selection Functions

Stellar information from the combined catalogue is subject to incompleteness arising from various factors, potentially introducing bias to the inferred density distribution. To mitigate this issue, it is crucial to correct the measured density distribution by accounting for the selection function of the stars. Following the approach outlined by Bovy (2017a), we estimate the selection function of the stars in the combined catalogue by comparing them to the assumed complete sample of stars from the reference catalogue (refer to subsection 2.1).

The selection function is defined as the probability of observing a star with specific properties in the catalogue or survey. Let $S(J,J-K_{s})$ represent the ratio between the number of stars in the combined catalogue and the number of stars in the reference catalogue with color $J-K_{s}$ and absolute magnitude $M_{J}$. The effective selection function is then defined as:

$$\mathfrak{S}(x,y,z)=\int\rho_{\textup{CMD}}(M_{J},c|x,y,z)S(M_{J},c)\textup{d}c\textup{d}M_{J},$$ (1)

where $c=J-K_{s}$, and $\rho_{\textup{CMD}}(M_{J},J-K_{s}|x,y,z)$ represents the density distribution of stars in the color-magnitude diagram (CMD) at spatial position $(x,y,z)$. Equation (1) can be seen as a conversion from color-magnitude space completeness to position space completeness. This equation is then carried to calculate the effective volume completeness, providing an estimate of the completeness for a given region $\Pi$, defined as:

$$\Xi(\Pi)=\frac{\int_{\Pi}\mathfrak{S}(x,y,z)\textup{d}^{3}x}{\int_{\Pi}\textup{d}^{3}x},$$ (2)

where in this case $\Pi$ represents a cylindrical volume with a radius of $150$ pc and varying height $z$ with some thickness $\Delta z$.

For this study, we make the assumption that the density distribution in the Color-Magnitude Diagram (CMD) is uniform. The values can thus be inferred from the Hertzsprung-Russell (HR) diagram of stars in the solar neighborhood, as illustrated in Fig. 2.

The effective volume completeness is employed to adjust the density distribution of stars inferred from the combined catalogue to the (supposedly) true density distribution of stars in the solar neighborhood. This correction involves dividing the density distribution of stars in the combined catalogue by the effective volume completeness, resulting in what is termed the corrected density distribution. Figure 3 displays samples from four groups of stars alongside their corresponding effective volume completeness. Notably, the effective volume completeness remains relatively constant for higher values of $|z|$, indicating that correction is essentially unnecessary for those regions.

Figure 3: The effective volume completeness function for stars in the solar neighborhood. This graph illustrates the behavior for four selected groups of stars, with other groups exhibiting similar trends.

3   MODELS

3.1 Distribution Function

In a self-gravitating system such as a galaxy, the distribution function (DF) of a star represents the probability density function of finding a star with a given position and velocity (or momentum), denoted as $f(\bm{x},\bm{v})$. In a collisionless system and under the condition of dynamical equilibrium, the DF adheres to the collisionless Boltzmann equation (CBE),

$$(\nabla_{\bm{x}}f)\cdot\bm{v}-(\nabla_{\bm{x}}\Phi)\cdot(\nabla_{\bm{v}}f) =0$$

$$(\nabla_{\bm{x}}f)\cdot\bm{v} =(\nabla_{\bm{x}}\Phi)\cdot(\nabla_{\bm{v}}f),$$ (3)

where $\Phi$ is the gravitational potential of the system in equilibrium. When considering the DF of a star in close proximity to the midplane, it can be approximated by separating the contributions into the vertical and horizontal directions,

$$f(\bm{x},\bm{v}) =Z(z,w)T(r,v_{r},\phi,v_{\phi}),$$ (4)

where $Z(z,w)$ is the DF in the vertical direction, and $T(r,v_{r},\phi,v_{\phi})$ is the DF in the horizontal direction. Equation (3) can be expressed as,

$$w\frac{\partial Z}{\partial z}-\frac{\partial\Phi}{\partial z}\frac{\partial Z}{\partial w} =0.$$ (5)

The differential equation (5) has a general solution $Z(z,w)=Z(E_{z})\equiv Z(w^{2}/2\Phi(z))$, where $E_{z}$ is the energy in the vertical direction. Thus, the DF in the vertical direction is a function of energy only. Assuming the star has a vertical velocity dispersion $f_{0}(w)$ in the midplane, the DF in the vertical direction can then be expressed as

$$Z(z,w) =Z(E_{z})$$

$$=Z(w^{2}/2+\Phi(z))$$

$$=Z(w_{0}^{2}/2+\Phi(z_{0}))$$

$$\propto f_{0}(w_{0})$$

$$\propto f_{0}\left(\sqrt{w^{2}+2[\Phi(z)-\Phi(z_{0})]}\right),$$ (6)

where, $z_{0}$ and $w_{0}$ represent the equivalent midplane position and velocity of the star that has the same energy as the star at $(z,w)$.

If we assume the vertical velocity dispersion $f_{0}$ follows a Gaussian distribution, then equation (6) can be solved analytically to obtain the DF of spatial vertical direction.

$$F_{z}(z) =\int_{-\infty}^{\infty}Z(z,w)\text{d}w$$

$$=A\exp\left(-\frac{\Phi(z)}{\sigma_{w}^{2}}\right),$$ (7)

$$\Rightarrow\nu(z) =\nu_{0}\exp\left(-\frac{\Phi(z)}{\sigma_{w}^{2}}\right),$$ (8)

where $A$ is a normalization constant, $\sigma_{w}$ is the vertical velocity dispersion, and $\Phi(z_{0})$ is chosen to be zero. If we only consider $N$ identical stars, then the vertical distribution (7) is proportional to the expected number density $\nu(z)$ of the stars, leading to (8).

The vertical velocity distribution can be calculated by integrating the DF in the spatial vertical direction,

$$F_{w}(w) =\int_{-\infty}^{\infty}Z(z,w)\text{d}z$$

$$=a^{*}f_{0}(w)$$

$$\Rightarrow\nu(w) =a\exp\left[-\frac{\left(w-w_{0}\right)^{2}}{2\sigma_{w}^{2}}\right],$$ (9)

where $a$ is a normalization constant. The vertical velocity distribution (9) is a Gaussian distribution with mean $w_{0}$ and standard deviation $\sigma_{w}$.

3.2 Gravitational Potential

The gravitational potential of a system can be calculated from the mass density distribution using Poisson’s equation. In a cylindrical volume around the Sun, Poisson’s equation can be expressed as

$$\nabla^{2}\Phi =\frac{\partial^{2}\Phi}{\partial z^{2}}+\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\Phi}{\partial r}\right)=4\pi G\rho$$

$$\frac{\partial^{2}\Phi}{\partial z^{2}} =4\pi G\rho_{\textup{eff}},$$ (10)

where the second term of equation (10) is referred to as the rotation curve term $\mathcal{R}$, which is effectively constant for a small volume around the Sun. The rotation curve term can be calculated as

$$\mathcal{R} =\frac{1}{4\pi G}\frac{1}{r}\frac{\partial}{\partial}\left(r\frac{\partial\Phi}{\partial r}\right)$$

$$=\frac{A^{2}-B^{2}}{2\pi G}$$ (11)

where $A$ and $B$ represent the Oort constants. The values of $A$ and $B$ are $15.3\pm 0.4$ km s${}^{-1}$ kpc${}^{-1}$ and $-11.9\pm 0.4$ km s${}^{-1}$ kpc${}^{-1}$, respectively (Bovy 2017b).

For the volume studied in this work, the rotation curve term is approximately constant and the gravitational potential $\Phi$ is only dependent on $z$. Therefore, the Poisson’s equation (10) can be written as a second-order ordinary differential equation with respect to $z$,

$$\frac{\textup{d}^{2}\Phi(z)}{\textup{d}z^{2}} =4\pi G\rho_{\textup{eff}}(\Phi,z),$$ (12)

where $\rho_{\textup{eff}}(\Phi,z)$ could be a function of $\Phi$ and $z$. Section 3.4 will discuss the mass model used in this work in detail.

3.3 Modified Poisson Equation

In Milgrom’s work 1983, he proposed a modification to Newtonian dynamics by introducing a new constant $a_{0}$, referred to as the acceleration parameter. This constant plays a crucial role in scenarios where an object experiences a low acceleration field. One approach to incorporating the acceleration parameter into Newtonian dynamics is by ensuring that the circular speed of an object in a low acceleration regime remains constant:

$$\text{Centripetal acc.}=\frac{V_{c}^{2}}{r}\propto\frac{1}{r}.$$ (13)

Such a modification can be achieved by adjusting the gravitating acceleration to follow

$$g=\sqrt{g_{N}a_{0}},$$ (14)

where $g_{N}$ represents the classical Newtonian gravitational acceleration.

Concerning this phenomenology, MOND introduces a different form of Poisson’s equation, as formulated by Bekenstein & Milgrom (1984) for AQUAL theory:

$$4\pi G\rho=\nabla\cdot\left[\mu\left(\frac{|\nabla\Phi|}{a_{0}}\right)\Phi\right],$$ (15)

where $a_{0}$ represents the MOND acceleration parameter constant, and $\mu$ is the MOND $\mu$-interpolating function. The $\mu$-interpolating function is defined as

$$\mu(x) =\begin{cases}1&\text{if }x\gg 1\\ x&\text{if }x\ll 1\end{cases}.$$ (16)

Equation (16) ensures that the modified Poisson’s equation reduces to the Newtonian Poisson’s equation when the acceleration is much larger than the MOND acceleration parameter $a_{0}$. In small volumes, where the total acceleration $|\nabla\Phi|$ is approximately constant, the modified Poisson’s equation (15) can be simplified to:

$$\nabla^{2}\Phi=\frac{4\pi G\rho_{\text{eff}}}{\mu_{0}},$$ (17)

where $\rho_{\text{eff}}$ is the effective mass density ($\rho-4\pi G\mathcal{R}$), and $\mu_{0}$ is the value of the interpolating function at the Sun’s position.

Both Poisson equations (12) and (17) are similar, with a minor difference represented by a factor of $\mu_{0}$. This value can be interpreted as the additional boost of the gravitational potential due to the MOND effect. In the dark matter hypothesis, this additional boost is generated by additional mass density from nonluminous dark matter. Essentially, the MOND effect can be mimicked by adding extra mass density to the system. The two equations can result in the exact same gravitational potential if the additional missing mass follows the baryonic mass distribution. Fortunately, most dark matter models predict that the dark matter distribution does not precisely follow the baryonic mass distribution (Navarro et al. 1997). As a result, the two equations should yield different gravitational potentials in theory. The extent of the difference between the two equations is determined by the value of $\mu_{0}$. If $\mu_{0}$ is close to 1, it becomes extremely challenging to distinguish between the hypotheses of dark matter and MOND.

The exact function of the $\mu$-interpolating function is not precisely known. However, for many studies focusing on the low acceleration regime, the specific form of the interpolating function is often not critical. In this study, where the acceleration is not low enough to be in the low acceleration regime, the exact function of the interpolating function becomes crucial for inferring the value of the acceleration parameter $a_{0}$. We consider two common interpolating functions: the simple and the standard. Table 3 provides the exact functions for these two interpolating functions.

Table 3: The formulas for the interpolating function $\mu(x)$.

Interpolating Function	$\mu(x)$
Simple	$\displaystyle\frac{x}{1+x}$
Standard	$\displaystyle\frac{x}{\sqrt{1+x^{2}}}$

3.4 Mass Model

The model employs in this study encompasses analytical functions for the vertical distribution of stars (8) and the vertical velocity distribution of stars (9). This model is derived by solving the Poisson equation (10 or 17) with the provided mass density distribution. The mass density distribution for this study is the sum of $N_{\textrm{b}}$ baryon components and a constant dark matter (DM) halo contribution in the midplane. The baryon mass density $\rho_{\textrm{b}}$ is described by the Bahcall model, which includes isothermal gas, stars, and stellar remnants (Bahcall 1984b, 1984c, 1984a).

$$\rho_{\textrm{b}}(z) =\sum_{i=1}^{N_{b}}\rho_{i}(0)\exp\left(-\frac{\Phi(z)}{\sigma_{z,i}^{2}}\right),$$ (18)

where $\rho_{i}(0)$ represents the midplane mass density of the $i$-th component at the midplane, and $\sigma_{z,i}$ denotes the vertical velocity dispersion of the $i$-th component. Table 4 provides the parameters of the Bahcall model utilized in this study.

Table 4: Baryonic mass components employed in this study. The table is extracted from table 2 of Buch et al. (2019).

Component	$\rho_{i,0}\left[\textrm{M}_{\odot}/\textup{pc}^{2}\right]$	$\sigma_{z,i}\left[\textup{km/s}\right]$
Molecular gas (H${}_{2}$)	$0.0104\pm 0.00312$	$3.7\pm 0.2$
Cold atomic gas (H$\,\text{\scriptsize I}(1)$)	$0.0277\pm 0.00554$	$7.1\pm 0.5$
Warm atomic gas (H$\,\text{\scriptsize I}(2)$)	$0.0073\pm 0.00070$	$22.1\pm 2.4$
Hot ionized gas (H II)	$0.0005\pm 0.00003$	$39.0\pm 4.0$
Giant stars	$0.0006\pm 0.00006$	$15.5\pm 1.6$
$M_{V}<3$	$0.0018\pm 0.00018$	$7.5\pm 2.0$
$3<M_{V}<4$	$0.0018\pm 0.00018$	$12.0\pm 2.4$
$4<M_{V}<5$	$0.0029\pm 0.00029$	$18.0\pm 1.8$
$5<M_{V}<8$	$0.0072\pm 0.00072$	$18.5\pm 1.9$
$M_{V}>8$ (M dwarfs)	$0.0216\pm 0.00280$	$18.5\pm 4.0$
White dwarfs	$0.0056\pm 0.00100$	$20.0\pm 5.0$
Brown dwarfs	$0.0015\pm 0.00050$	$20.0\pm 5.0$

The DM halo contribution is assumed to be constant in the vicinity of the Sun, so the total mass density distribution is given by

$$\rho(z) =\rho_{\textrm{b}}(z)+\rho_{\textrm{\scriptsize DM}}.$$ (19)

Based on the mass density distribution, the gravitational potential is calculated by solving the Poisson equation (10) or (17). This potential is then used to compute the vertical distribution of stars (8) and the vertical velocity distribution of stars (9).

Solving the Poisson equation with this mass model is not trivial because it is not analytically solvable or separable. Numerical methods are needed to solve the Poisson equation. Fundamentally, from the Poisson equation, given mass model $\rho$, there is a unique gravitational potential $\Phi$.

In this study, we explore three hypotheses: the baseline hypothesis involves Classic Newtonian gravity without Dark Matter (NO), while the other two hypotheses consider Newtonian gravity with dark matter (DM) and Modified Newtonian Dynamics (MOND). The DM model incorporates the DM halo contribution into the mass density distribution, while the MOND model relies solely on baryonic mass contributions. Despite the exclusion of the DM halo contribution in the MOND model’s equation (19), it introduces the MOND effect into the Poisson equation (17). Equation (17) can be perceived as the traditional Poisson equation with baryonic mass density, scaled by a factor of $1/\mu_{0}$ times the baryonic matter. Essentially, for our specific study, the MOND model is equivalent to a scenario without a DM model using classical Newtonian gravity, albeit with adjusted mass density.

4   METHODS

4.1 Density Inference

The most straightforward method for inferring the density distribution of stars from their spatial distribution is to count the number of stars within a specified volume. While this approach is rudimentary, its simplicity is a powerful tool, particularly when dealing with a large number of stars. However, a drawback of this method is the omission of uncertainties in star positions, which may be significant in certain cases.

Previous attempts, as demonstrated by Buch et al. (2019); Xia et al. (2016), Sivertsson (2018), among others, have aimed to infer density uncertainty by employing Poisson error or simple standard deviation within a given volume. While this approach is superior to having no uncertainty estimation, it often assumes that the distribution follows a Gaussian (or even log-Gaussian) shape, which is not always accurate. For a large number of stars, the Poisson distribution can be well approximated by a Gaussian distribution. However, when dealing with a small number of stars, assuming a Gaussian distribution might lead to inaccuracies. In such cases, the distribution should be considered back to Poisson distribution.

While the Poisson error in a given volume accounts for uncertainty, it does not consider the individual star’s position uncertainty. For datasets with high-quality parallax data, individual star position uncertainties may be negligible. However, in our study, we aim to develop a method applicable to any type of data, including those with poor quality and significant uncertainty. Appendix A outlines the methodology for inferring the density distribution of objects from their ‘position’. Using this approach, we can infer the density distribution of stars from their spatial distribution, taking into account the uncertainty in their positions. In other words, for a given distance $z$ from the midplane, there is a probability density function $P(\nu,z)$ representing the likelihood of the density being $\nu$.

We employ the technique proposed by Bailer-Jones (2015) to obtain the probability distribution function of the distance $r$ with a constant volume density prior:

$$p(r|\varpi,\sigma_{\varpi})\propto r^{2}\exp\left[-\frac{1}{2\sigma_{\varpi}^{2}}\left(\varpi-\frac{1}{r}\right)^{2}\right].$$ (20)

Subsequently, we convert this distribution to the vertical distance $z$,

$$p(z|\varpi,\sigma_{\varpi},b)\propto z^{2}\exp\left[-\frac{1}{2\sigma_{\varpi}^{2}}\left(\varpi-\frac{\sin b}{z}\right)^{2}\right],$$ (21)

where $(\varpi,\sigma_{\varpi})$ represent the parallax and its error, $z$ denotes the vertical distance from the Galactic midplane, and $b$ is the Galactic latitude.

A similar approach is used to determine the probability distribution function of vertical velocity. The formula used to calculate $w$ from its proper motion and radial velocity is given by:

$$w(\mu_{b},b,v_{r},r)=\kappa\mu_{b}r\cos b+v_{r}\sin b,$$ (22)

where $\kappa=4.74$ km s${}^{-1}$ pc${}^{-1}$ arcsec yr${}^{-1}$, $\mu_{b}$ is the proper motion in the galactic latitude direction, and $v_{r}$ is the radial velocity.

We generate 100,000 samples of ($\varpi,\mu_{b},v_{r}$) for each star, then calculate their $z$ and $w$ values. Subsequently, we bin both $z$ and $w$ for all star groups and calculate their weights $W$ according to the methods outlined in Appendix A. Fig 4 shows one such inference for the tracer of the first group (A0-A8).

4.2 Likelihood Function

The likelihood function is divided into two parts: the likelihood function of the vertical number density $z$ and the likelihood function of the vertical velocity dispersion $w$. Given a set of parameters $\zeta$ that describe the vertical distribution of stars, the likelihood function of finding $N(z_{i}|\zeta)$ stars in the $i$-th bin is expressed as

$$\mathcal{L}(N(z_{i}|\zeta)) =\sum_{n=0|z=z_{i}}^{\infty}W(n)\textup{Poisson}\left(n|N(z_{i}|\zeta)\right),$$ (23)

where $W(n)$ represents the weight of finding $n$ stars in the $i$-th bin, calculated from the data collected in Section 2. The value of $N(z_{i}|\zeta)$ is calculated using (8). The term $\textup{Poisson}(n|N)$ is the probability density function of the Poisson distribution, representing the likelihood of finding $n$ stars given the expected number of stars $N$. For a more detailed derivation, refer to Appendix A.

The same function is used to calculate the likelihood function of finding $N(w_{i}|\xi)$ stars in the $i$-th bin of the velocity distribution, where $\xi$ is a set of parameters defining the vertical velocity distribution of the stars. Therefore,

$$\mathcal{L}(N(w_{i}|\xi)) =\sum_{n=0|z=z_{i}}^{\infty}W(n)\textup{Poisson}\left(n|N(w_{i}|\xi)\right).$$ (24)

The total likelihood function is the product of the likelihood function for the vertical number density and the likelihood function for the vertical velocity dispersion. The complete set of free parameters $\theta$ encompasses both the combined dynamics parameters $\zeta$ and the kinematics parameters $\xi$.

4.3 Priors

The prior distribution of the parameters $\zeta$ and $\xi$ is summarized in Table 6 and Table 5, respectively. We assume that for each group, there are actually two populations of stars: one with a broader velocity dispersion and another with a narrower velocity dispersion. Most likely, these two populations correspond to the thin disc and the thick disc and/or halo. The thin disc is assumed to be the dominant population, so its contribution is assumed to be higher than that of the contaminant component. This assumption is based on the fact that the thin disc is the dominant population in the solar neighborhood. Adopting this approach, we observe that the fitting curve is improved compared to using only one population.

Only two types of prior distributions are used: uniform distribution and normal distribution. In Table 6 and Table 5, the “Constraints” column is defined as lower and upper bounds $[a,b]$ if the distribution is uniform, and mean and standard deviation $[\mu,\sigma]$ if the distribution is normal.

We also impose restrictions on the joint prior distribution of the parameters $\sigma_{w,1}$ and $\sigma_{w,2}$ to be $\sigma_{w,1}<\sigma_{w,2}$ and on the parameters $a_{1}$ and $a_{2}$ to be $a_{1}>a_{2}$. Parameter $a$ is the sum of $a_{1}$ and $a_{2}$. These restrictions are implemented based on the assumption that the thin disc is the dominant population and that the thin disc has a narrower velocity dispersion than the contaminant population.

Table 5: Summary of kinematic priors, with units indicated in parentheses corresponding to the argument of logarithm.

Parameter	Distribution	Contraints	Unit
$w_{0}$	Uniform	$[-15,0]$	km/s
$\log\sigma_{w,i}$	Uniform	$[\ln 3,\ln 50]$	[km/s]
$\log a$	Normal	$[\ln\textup{max}(N(w_{i})),2]$	-
$\log a_{2}$	Uniform	$[\ln\textup{max}(N(w_{i}))-5,\ln\textup{max}(N(w_{i}))]$	-

Table 6: Summary of dynamic priors. The constraints for baryon parameters are referenced to Table 4.

Parameter	Distribution	Contraints	Unit
$\rho_{\textrm{b},i}$	Normal	Table 4	M${}_{\odot}$ pc${}^{-3}$
$\sigma_{z,i}$	Normal	Table 4	km s${}^{-1}$
$\rho_{\textup{\tiny DM}}$	Uniform	[-0.05, 0.15]	M${}_{\odot}$ pc${}^{-3}$
$\ln\mu_{0}$	Uniform	[$\ln 0.1,\ln 2.0$]	-
$\ln\nu_{0}$	Normal	[$\ln\textup{max}(\vec{N}(z_{i}))$, 3]	-
$\mathcal{R}$	Normal	[3.4, 0.6]$\times 10^{-3}$	M${}_{\odot}$ pc${}^{-3}$
$z_{0}$	Uniform	[-150, 150]	pc

4.4 MCMC

We make use of the Markov Chain Monte Carlo (MCMC) method to sample the posterior distribution of the parameters $\theta$. The MCMC simulations are conducted using hammer-and-sample by Adam Reichold, a Rust implementation of the affine invariant ensemble sampler introduced by Goodman & Weare (2010), which is a reimplementation of emcee Python library by Foreman-Mackey et al. (2013).

For each group, the MCMC is run for 50,000 steps with 10 $\times$ ndim walkers, where ndim $=33$ for DM or MOND model, and ndim $=32$ for classic NO model. The first 4,000 iterations are discarded as burn-in, and the samples are thinned by a factor of 20.

The MCMC program is implemented in the Rust programming language to enhance performance and reliability. With 33 free parameters for each group in the DM model, spanning 15 groups based on infrared color, and considering 2 gravitational models, a substantial amount of computation is involved. Although the MCMC program is easy to implement, optimizing it can be quite challenging. Rust demonstrates superior performance, outperforming even well-written Python code. While Rust may not be as user-friendly as Python, and its ecosystem is not as mature in Science and Engineering, the significant performance gains justify the effort. Despite this, the program is also ported to Python for the sake of easier data analysis and visualization.

4.5 Solving the Poisson Equation

The solution to the Poisson equation is achieved using the Dormand-Prince method of order 5(4) with dense output of order 4 (Dopri5). This method is implemented in the Rust crate ode_solvers by Sylvain Renevey. The choice of this method is driven by its speed, ease of use, and direct integration into Rust. The output from the integration seamlessly interfaces with the rest of the program in Rust, eliminating the need for data structure conversion or cross-language interfaces. The native code compilation enhances the program’s performance blazingly fast.

The general idea of the method involves providing the initial conditions $\Phi(0)=0$ and $\partial_{z}\Phi(0)=0$, along with a set of parameters $\theta$. The method then performs the integration of the Poisson equation at discrete values of $z$, yielding the values of $\Phi(z)$ and $\partial_{z}\Phi(z)$ at those points. The resulting integration is subsequently interpolated using a simple linear interpolation method to obtain the value of $\Phi(z)$ at any arbitrary $z$. This calculated $\Phi(z)$ value is then used to determine the vertical number density of stars using (8).

4.6 Bayes Factor Estimation

The computation of the Bayes factor from the MCMC chain involves employing the approach presented in Robert & Wraith (2009). The harmonic means method is utilized to estimate the evidence term for each MCMC simulation, given by

$$E\approx N\left(\sum_{i=1}^{N}\frac{\phi(\theta_{i})}{\mathcal{L}(\theta_{i})p_{0}(\theta_{i})}\right)^{-1},$$ (25)

where $N$ denotes the number of samples obtained from the MCMC simulations, and $\phi(\theta_{i})$ is any arbitrary probability distribution function over $\theta$. While a simple harmonic means method estimates the evidence term by selecting $\phi=p_{0}$, it can be risky due to the instability of the approximation (Newton & Raftery 1994). In our case, we implement the uniform distribution $\phi(\theta)$ with boundaries around the peak of the posterior distribution $p(\theta)$, following the strategy outlined in Robert & Wraith (2009), which stabilizes the approximation. It is crucial to note that for this method to be effective, the proper prior distribution function $p_{0}$ must be used. The Bayes factor is defined as the ratio between the evidence terms of different models, $B_{12}=E_{1}/E_{2}$. We calculate the approximation 10 times, then take the mean and the standard deviation as the estimated value.

5   RESULTS AND DISCUSSION

After running the MCMC for each group and model, the results are in the form of a 3-dimensional array with respect to the number of walkers, number of parameters, and number of steps. These results are then flattened (combining all walkers) to obtain the full samples of the posterior distribution functions $\theta$. The dark matter density $\rho_{\textrm{\tiny DM}}$ and boost interpolating function value $\mu_{0}$ are extracted by constructing their histograms. Appropriate functions are used to obtain their central values and credible intervals.

5.1 Dark Matter Density

The full set of estimated parameters $\theta$ is best-fitted to the inferred density profile using MCMC (Fig. 4). This reveals that the contribution of dark matter in the halo is approximately on the order of $0.01\,\textrm{M}_{\odot},\textrm{pc}^{-3}$. The results exhibit fluctuations and are not single-valued, as illustrated in Fig. 5. The figure also displays the inferred total baryonic matter density profile, which aligns with the initial values provided in Table 4. This indicates that the fitting process accurately captures the contribution of the known baryonic matter.

Figure 4: Density inference (error-bars) and best-fitting curve (shaded area) to the vertical profile density and vertical velocity dispersion for one of the tracer groups (A0-A8). All uncertainties are presented in 68, 90, and 95 percent credible intervals.

In broad terms, the inferred dark matter densities from A-type and G-type stars exhibit agreement, hovering around $\sim 0.01-0.03$ M${}_{\odot}$ pc${}^{-3}$, while the value for F-type stars is notably higher, reaching up to $\sim 0.07$ M${}_{\odot}$ pc${}^{-3}$. The combined dark matter density for all groups is $(3.71^{+0.39}_{-0.36})\times 10^{-2}$ M${}_{\odot}$ pc${}^{-3}$, while if F-type stars are excluded, the combined infered dark matter density becomes $(1.83^{+0.53}_{-0.52})\times 10^{-2}$ M${}_{\odot}$ pc${}^{-3}$. Theoretically, the dark matter density should be consistent across all spectral types, as the dark matter gravitational potential contribution is expected to be independent of the spectral type of stars. The implication of this result could stem from either the unsuitability of the dark matter model employed in this study or the violation of some assumptions made in the analysis. Let us discuss the possible causes of this result.

Firstly, it’s important to note that the dark matter model employed in this study is relatively simple, assuming a constant density up to 200 pc above and below the galactic plane. However, simulations of Milky Way-like galaxies, as indicated by Garbari et al. (2011), suggest that the dark matter density remains approximately constant up to 1 kpc above and below the galactic plane, even after 4 Gyrs of evolution. While the simulation may not perfectly reflect the real Galaxy, it provides valuable insights into the nature of dark matter distribution. Consequently, it seems unlikely that refining the dark matter density model would significantly alter the obtained results.

Secondly, among the various assumptions made in this study, two stand out as potential contributors to the observed peculiar results. The first assumption involves the baryon mass distribution, following the Bahcall model, which consists of several isothermal components. The isothermality assumption implies that the velocity dispersion remains constant up to 200 pc. Fortunately, simulations of Milky Way-like galaxies, as suggested by Garbari et al. (2011), indicate that the velocity dispersion is approximately constant up to 1 kpc above and below the galactic plane. Thus, this assumption is likely valid. The second assumption that raises suspicion is the assumption of dynamical equilibrium, implying that the stars under consideration maintain the same distribution function after years of interaction and mixing. Achieving dynamical equilibrium requires that the stars are old enough to interact and mix but not so old that they have been evaporated from the galaxy. Furthermore, dynamical equilibrium can only be sustained if there are no significant perturbations to the system, such as a merger with other galaxies (Malhan et al. 2022). This might be the primary cause of the fluctuating results in the inferred dark matter density for different groups of stars within distinct color ranges.

The lack of convergence in the results for dark matter density strongly suggests that the sampled data is not in dynamical equilibrium. This finding aligns with the results reported by Buch et al. (2019), indicating that the dark matter density differs for various spectral types of stars. To address this issue, a modification to the collisionless Boltzmann equation may be necessary to account for the effects of disequilibrium. One potential approach involves introducing external wave-like perturbations to the distribution function, as proposed by (Banik et al. 2017). This presents an intriguing avenue for future research.

Figure 5: Inferred dark matter density for each group of stars. The vertical dashed line separates different spectral types, while the horizontal dashed line represents the initial total baryonic matter density from Table 4, with the gray diagonal shaded area indicating the uncertainty. The red error bar plot represents the inferred baryonic matter density, and the blue error bar plot depicts the inferred dark matter density from the fitting process. All uncertainties are presented in 68, 90, and 95 per cent credible intervals.

5.2 Acceleration Parameters

Fitting the data without including dark matter contribution reveals that the inferred baryonic matter density consistently exceeds the initial value from Table 4. This suggests the presence of another matter component, such as dark matter, or the possibility that the additional matter is an effect of modified gravity.

Equation (17) shows the relationship between baryonic matter and the additional $\mu_{0}$ boost from MOND theory. Our analysis yields a value of $\mu_{0}$ around $\sim 0.64\pm 0.03$ for the Sun’s location in the galaxy, without delving into the specifics of MOND theory.

A similar method to test MOND theory was conducted by Lisanti et al. (2019), revealing that the value of the $\nu$-interpolating function is $1.44^{+0.13}_{-0.08}$. Our equivalent $\nu$-interpolating function value is $1.56\pm 0.06$, which falls within the uncertainty range. This suggests that our analysis results are consistent with their findings.

MOND theory places particular emphasis on determining the value of the acceleration parameter $a_{0}$. Converting the value of $\mu_{0}$ into acceleration parameters $a_{0}$ requires a specific interpolating function. We use two interpolating functions: the simple one and the standard one, as outlined in Table 3. The resulting value of $a_{0}$ for the simple interpolating function is approximately $(1.26^{+0.14}_{-0.13})\times 10^{-10}$ m s${}^{-2}$, while for the standard interpolating function, it is around $(2.69^{+0.19}_{-0.17})\times 10^{-10}$m s${}^{-2}$. Fig. 6 shows the probability density function for simple and standard interpolating function. An intriguing observation is that the values of $a_{0}$ for these two interpolating functions differ significantly. Upon closer inspection, it is noted that the uncertainties of both results do not overlap, suggesting the potential to distinguish between the two models.

Previous studies by Gentile et al. (2011) and Begeman et al. (1991) found that the value of $a_{0}$ is approximately $1.2\times 10^{-10}$ m s${}^{-2}$ by fitting the rotation curve of galaxies. Using the methodology of our approach, constraining the gravitational potential by proxy of the local distribution function, we found a similar value. If MOND is valid, the simple interpolating function is more likely to be true to achieve the same result from a global galactic-scale approach.

Figure 6: The acceleration parameter probability density derived from the infered $\mu_{0}$. The gray curve in inset plot are the infered $\mu_{0}$ for each groups and dark bold curve is the combined values for all groups, around $0.64$.

Adopting the standard interpolating function, the result of $a_{0}\sim 1.2\times 10^{-10}$ m s${}^{-2}$ is consistent with other studies employing MOND on Galactic scales. This finding indicates that even in the solar neighborhood, the prospect to observe the MOND effect is significant and compelling. Further comprehensive investigations on this scale could shed light on the true nature of MOND and Dark Matter.

5.3 Comparisons

Bayes factors (BF) serve as a tool for comparing two or more models based on the observed data. The Bayes factor represents the integrated posterior probability of a model across all potential parameter values, often referred to as the “evidence” term in Bayesian statistics. A higher Bayes factor indicates that a model is more suitable for explaining the observed data.

The three hypotheses being compared are $M_{\textrm{no}}=$ (Newtonian Dynamics) with no dark matter, $M_{\textrm{dm}}=$ Dark Matter, and $M_{\textrm{mond}}=$ MOND. Fig. 7 presents the complete results. The mean $\log$ BF values with respect to $M_{\textrm{no}}=$ for all groups of stars are $1.80$ (dm-no) and $1.70$ (mond-no), suggesting that any additional matter from the baryonic dark matter hypothesis (no) is strongly unlikely. This trend is more pronounced for stars of type F, where the value spikes to $6$. Apart from these spikes, the $\log$ BF values are close to zero, indicating that the three hypotheses are equally likely to be true. The only strong conclusion that can be drawn is Newtonian dynamics without including non-baryonic dark matter is not feasible.

Compare to the Bayes factor between the Dark Matter and MOND hypotheses; the two are neck and neck. The $\log$ BF fluctuates around zero, with a mean of only 0.10 in favor of dark matter. This result suggests that both hypotheses are equally plausible in explaining the data we have. From our analysis, we cannot definitively say which hypothesis is more favored.

One way to better discriminate between the two hypotheses is by incorporating more data from stars at high altitudes ($|z|>200$ pc). This approach aims to reveal potential differences in the predictions of the distribution function between the two hypotheses. However, a challenge arises for large altitudes as the assumption that the distribution function is separable becomes invalid. Consequently, solving the full form of the distribution function is necessary (see (3)). This poses a complex and challenging problem to address.

Other studies, such as the one conducted by Lisanti et al. (2019), have employed statistical methods to compare the Dark Matter and MOND hypotheses, specifically by calculating the Bayesian Information Criterion (BIC). Their findings strongly favor the Dark Matter hypothesis based on the difference in BIC values ($\Delta\textup{BIC}$), which contrasts with our results. In our analysis, the Bayes factors for the Dark Matter hypothesis are only marginally higher than those for the MOND hypothesis. However, it is important to note that our results do not suggest strong evidence against the MOND hypothesis, with mean $\log$ BF values of $0.1$. According to Kass & Raftery (1995), a log Bayes factor of less than 0.5 is not even worth more than a bare mention. More data and robust methods are needed to provide conclusive answers.

Figure 7: Bayes factors for each group of stars for the three hypotheses. (Top) Bayes factors of Dark Matter and MOND hypotheses with respect to NO hypothesis. The bayes factor for MOND-NO points are slightly shifted to the left for clarity.(Bottom) Bayes factors comparing the Dark Matter hypothesis to the MOND hypothesis.

6   CONCLUSION

In this study, we have developed a method to estimate the dark matter density in the solar neighbourhood by leveraging the vertical distribution of stars obtained from Gaia, RAVE, APOGEE, GALAH, and LAMOST catalogues. The data collected from these various catalogues are dominated by Gaia DR3. In principle, the analysis done in this study is unaffected by the addition of other catalogues.

Our approach utilises both the vertical height distribution and vertical velocity distribution of stars as the vertical distribution function. In turn, the distribution function is directly correlated with the Galactic gravitational potential through CBE, which dictates the mass model for our hypotheses. By proxy of the distribution function, the dark matter density and the acceleration parameter $a_{0}$ is inferred.

Our analysis reveals a dark matter density in the solar neighbourhood of approximately $(0.01\textrm{--}0.03)\textrm{ M}_{\odot}$ pc${}^{-3}$, and an acceleration parameter $a_{0}$ of $1.26^{+0.14}_{-0.13}\times 10^{-10}$ m s${}^{-2}$ for a simple interpolating function. These fluctuations could be the results of disequilibria of tracers being used, suggesting challenges to improve the approach for more robust results in future studies.

Bayesian analysis, through the calculation of Bayes factors, presents strong evidence against classical Newtonian gravity without additional missing mass. However, when comparing the two competing hypotheses in explaining the missing mass, we found a marginally higher likelihood for the dark matter hypothesis compared to the Modified Newtonian Dynamics (MOND) hypothesis, with $\log\textrm{BF}=0.1$. It is important to note, however, that the Bayes factors provide no substantial support for the dark matter hypothesis, falling short of conclusive evidence to reject the MOND hypothesis outright.

Data Availability

The original raw data are available from publicly availabe source at their respective repositories. The reduced data underlying this article will be shared upon reasonable request to the corresponding author.

Appendix A Density Inference

Let $x$ be an observed value with some uncertainty sampled from a distribution $\mathcal{T}$, with mean $\mu$ and standard deviation $\sigma$. The mean value $\mu$ can be inferred using Bayes’ theorem:

$$p(\mu|x,\sigma) =\frac{p(x|\mu,\sigma)p_{0}(\mu|\sigma)}{E(x,\sigma)}$$

$$\propto\mathcal{T}(x|\mu,\sigma),$$ (26)

where $p_{0}(\mu|\sigma)\propto 1$ is assumed. If there is a set of observed values $\{\bm{x},\bm{\sigma}\}$ that are independent and identically distributed (iid), the probability density function of the mean value $\bm{\mu}$ is

$$p(\bm{\mu}|\bm{x},\bm{\sigma}) \propto\prod_{i=1}^{N}\mathcal{T}(x_{i}|\mu_{i},\sigma_{i}).$$ (27)

Consider an interval $x$ and $x+\Delta x$. Due to the noisy nature of the data, the number of observed values in that interval is a random variable. Let $N$ be the number of data points with mean $\mu$ value within the interval. The probability density function of $N$ is given by

$$p(N|\bm{x},\bm{\sigma}) =\int p(N,\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$$

$$=\int p(N|\bm{\mu},\bm{x},\bm{\sigma})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$$

$$=\int p(N|\bm{\mu})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$$

$$=\int\mathbb{I}(\bm{\mu})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu},$$ (28)

where $\mathbb{I}(\bm{\mu})$ is one when the combination of $\bm{\mu}$ results in $N$ observed values in the interval, and zero otherwise. In reality, the random variable $\mu$ is distributed according to a Poisson distribution. Let $\Lambda$ be the expected number of observed values of $\mu$ in the interval; then, the probability density function of $\Lambda$ is given by

$$p(\Lambda|\bm{x},\bm{\sigma}) =\sum p(\Lambda,N_{i}|\bm{x},\bm{\sigma})$$

$$=\sum p(\Lambda|N_{i},\bm{x},\bm{\sigma})p(N_{i}|\bm{x},\bm{\sigma})$$

$$=\sum p(\Lambda|N_{i})\int\mathbb{I}(\bm{\mu})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$$

$$=\sum W(N_{i})p(\Lambda|N_{i}),$$ (29)

where $p(\Lambda|N_{i})\propto p(N_{i}|\Lambda)$ is the probability density function of the Poisson distribution with mean $\Lambda$ given an uninformative prior, and $W(N_{i})$ is the weight of finding $N_{i}$ observed values in the interval. The weight $W(N_{i})$ is extremely hard to calculate from the integral in equation (28), so in practice we use the following approximation:

$$W(N_{i})\approx\frac{1}{S}\sum_{j=1}^{S}\mathbb{I}(\bm{\mu}_{j}),$$ (30)

where we sample $\bm{\mu}$ from $p(\bm{\mu}|\bm{x},\bm{\sigma})$ $S$ times.

Equation (29) is computationally inefficient due to the large number of possible values for $N_{i}$. To address this, we defined a new distribution function called the modified Poisson distribution to approximate the probability density function of $\Lambda$:

$$\textrm{Poisson}^{*}(N|\Lambda)=\frac{\Lambda^{N}e^{-\Lambda}}{\Gamma(N+1)},$$ (31)

where $\Gamma(x)$ is the gamma function. Our tests show that the modified Poisson distribution is a good approximation for equation (29). Thus, we express the probability density function of $\Lambda$ as

$$p(\Lambda|\bm{x},\bm{\sigma}) \approx\textrm{Poisson}^{*}(N_{\textrm{eff}}|\Lambda),$$ (32)

where $N_{\textrm{eff}}$ represents the effective number of observed values in the interval. This number is not necessarily an integer. This choice minimizes the mean square error between equation (29) and equation (32), ensuring both computational efficiency and a close approximation to the true value.