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

M. A. Syaifudin,1 M. I. Arifyanto,2,3 H. R. T. Wulandari2,3, F. A. M. Mulki2,3
1Astronomy Master Program, Bandung Institute of Technology, Jl. Ganesa No. 10,
40135, Indonesia,
2Department of Astronomy, Bandung Institute of Technology, Jl. Ganesa No. 10,
40135, Indonesia,
3Bosscha Observatory, Bandung Institute of Technology, Jl. Peneropongan
Bintang No. 1, 40391, Indonesia

###### 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⊙ pc-3, although some studies show higher values, reaching
up to $\sim 0.1$ M⊙ 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
archive111https://gea.esac.esa.int/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)222We use version 2022.04.16, https://www.pas.rochester.edu/
emamajek/EEM_dwarf_UBVIJHK_colours_Teff.txt.

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:

$\displaystyle\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:

$\displaystyle\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),

$\displaystyle(\nabla_{\bm{x}}f)\cdot\bm{v}-(\nabla_{\bm{x}}\Phi)\cdot(\nabla_{\bm{v}}f)$
$\displaystyle=0$ $\displaystyle(\nabla_{\bm{x}}f)\cdot\bm{v}$
$\displaystyle=(\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,

$\displaystyle f(\bm{x},\bm{v})$
$\displaystyle=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,

$\displaystyle w\frac{\partial Z}{\partial z}-\frac{\partial\Phi}{\partial
z}\frac{\partial Z}{\partial w}$ $\displaystyle=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

$\displaystyle Z(z,w)$ $\displaystyle=Z(E_{z})$
$\displaystyle=Z(w^{2}/2+\Phi(z))$ $\displaystyle=Z(w_{0}^{2}/2+\Phi(z_{0}))$
$\displaystyle\propto f_{0}(w_{0})$ $\displaystyle\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.

$\displaystyle F_{z}(z)$
$\displaystyle=\int_{-\infty}^{\infty}Z(z,w)\text{d}w$
$\displaystyle=A\exp\left(-\frac{\Phi(z)}{\sigma_{w}^{2}}\right),$ (7)
$\displaystyle\Rightarrow\nu(z)$
$\displaystyle=\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,

$\displaystyle F_{w}(w)$
$\displaystyle=\int_{-\infty}^{\infty}Z(z,w)\text{d}z$
$\displaystyle=a^{*}f_{0}(w)$ $\displaystyle\Rightarrow\nu(w)$
$\displaystyle=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

$\displaystyle\nabla^{2}\Phi$ $\displaystyle=\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$
$\displaystyle\frac{\partial^{2}\Phi}{\partial z^{2}}$ $\displaystyle=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

$\displaystyle\mathcal{R}$ $\displaystyle=\frac{1}{4\pi
G}\frac{1}{r}\frac{\partial}{\partial}\left(r\frac{\partial\Phi}{\partial
r}\right)$ $\displaystyle=\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$,

$\displaystyle\frac{\textup{d}^{2}\Phi(z)}{\textup{d}z^{2}}$
$\displaystyle=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:

$\displaystyle\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

$\displaystyle 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:

$\displaystyle 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

$\displaystyle\mu(x)$ $\displaystyle=\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:

$\displaystyle\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).

$\displaystyle\rho_{\textrm{b}}(z)$
$\displaystyle=\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 (H2) | $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

$\displaystyle\rho(z)$
$\displaystyle=\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:

$\displaystyle 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$,

$\displaystyle 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:

$\displaystyle 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

$\displaystyle\mathcal{L}(N(z_{i}|\zeta))$
$\displaystyle=\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,

$\displaystyle\mathcal{L}(N(w_{i}|\xi))$
$\displaystyle=\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⊙ pc-3
$\sigma_{z,i}$ | Normal | Table 4 | km s-1
$\rho_{\textup{\tiny DM}}$ | Uniform | [-0.05, 0.15] | M⊙ 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⊙ 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-sample333https://crates.io/crates/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_solvers444https://crates.io/crates/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

$\displaystyle 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⊙ pc-3, while the
value for F-type stars is notably higher, reaching up to $\sim 0.07$ M⊙ pc-3.
The combined dark matter density for all groups is
$(3.71^{+0.39}_{-0.36})\times 10^{-2}$ M⊙ 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⊙ 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.

## Acknowledgements

We thanked the Faculty of Mathematics and Natural Sciences of Bandung
Institute of Technology for financial support via the PPMI 2022 funding
program.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

This publication makes use of data products from the Two Micron All Sky
Survey, which is a joint project of the University of Massachusetts and the
Infrared Processing and Analysis Center/California Institute of Technology,
funded by the National Aeronautics and Space Administration and the National
Science Foundation.

Funding for Rave has been provided by: the Leibniz Institute for Astrophysics
Potsdam (AIP); the Australian Astronomical Observatory; the Australian
National University; the Australian Research Council; the French National
Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the
European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale
di Astrofisica at Padova; the Johns Hopkins University; the National Science
Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie
University; the Netherlands Research School for Astronomy; the Natural
Sciences and Engineering Research Council of Canada; the Slovenian Research
Agency; the Swiss National Science Foundation; the Science & Technology
Facilities Council of the UK; Opticon; Strasbourg Observatory; the
Universities of Basel, Groningen, Heidelberg and Sydney. This work has made
use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope LAMOST) is a National Major Scientific Project built by the Chinese
Academy of Sciences. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences.

This work made use of the Third Data Release of the GALAH Survey (Buder et al.
2021). The GALAH Survey is based on data acquired through the Australian
Astronomical Observatory, under programs: A/2013B/13 (The GALAH pilot survey);
A/2014A/25, A/2015A/19, A2017A/18 (The GALAH survey phase 1); A2018A/18 (Open
clusters with HERMES); A2019A/1 (Hierarchical star formation in Ori OB1);
A2019A/15 (The GALAH survey phase 2); A/2015B/19, A/2016A/22, A/2016B/10,
A/2017B/16, A/2018B/15 (The HERMES-TESS program); and A/2015A/3, A/2015B/1,
A/2015B/19, A/2016A/22, A/2016B/12, A/2017A/14 (The HERMES K2-follow-up
program). We acknowledge the traditional owners of the land on which the AAT
stands, the Gamilaraay people, and pay our respects to elders past and
present. This paper includes data that has been provided by AAO Data Central
(datacentral.org.au).

Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P.
Sloan Foundation, the Heising-Simons Foundation, the National Science
Foundation, and the Participating Institutions. SDSS acknowledges support and
resources from the Center for High-Performance Computing at the University of
Utah. The SDSS web site is www.sdss.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration, including the Carnegie Institution for
Science, Chilean National Time Allocation Committee (CNTAC) ratified
researchers, the Gotham Participation Group, Harvard University, Heidelberg
University, The Johns Hopkins University, L’Ecole polytechnique fédérale de
Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-
Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für
Extraterrestrische Physik (MPE), Nanjing University, National Astronomical
Observatories of China (NAOC), New Mexico State University, The Ohio State
University, Pennsylvania State University, Smithsonian Astrophysical
Observatory, Space Telescope Science Institute (STScI), the Stellar
Astrophysics Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Illinois
at Urbana-Champaign, University of Toronto, University of Utah, University of
Virginia, Yale University, and Yunnan University.

## 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.

## References

  * Rubin et al. [1980] V. C. Rubin, Jr. Ford, W. K., and N. Thonnard.  Rotational properties of 21 SC galaxies with a large range of luminosities and radii, from NGC 4605 (R=4kpc) to UGC 2885 (R=122kpc).  _ApJ_ , 238:471–487, June 1980.  doi: 10.1086/158003.
  * Read [2014] J. I. Read.  The local dark matter density.  _J. Phys. G: Nucl. Phys_ , 41(6):063101, June 2014.  doi: 10.1088/0954-3899/41/6/063101.
  * Milgrom [2001] Mordehai Milgrom.  The shape of ‘dark matter’ haloes of disc galaxies according to mond.  _MNRAS_ , 326(4):1261–1264, 2001.  doi: 10.1111/j.1365-2966.2001.04653.x.
  * Nipoti et al. [2007] Carlo Nipoti, Pasquale Londrillo, HongSheng Zhao, and Luca Ciotti.  Vertical dynamics of disc galaxies in modified Newtonian dynamics.  _MNRAS_ , 379(2):597–604, August 2007.  doi: 10.1111/j.1365-2966.2007.11835.x.
  * Oort [1932] J. H. Oort.  The force exerted by the stellar system in the direction perpendicular to the galactic plane and some related problems.  _Bull. Astron. Inst. Neth._ , 6:249, August 1932.
  * Zwicky [1933] F. Zwicky.  Die Rotverschiebung von extragalaktischen Nebeln.  _Helv. Phys. Acta_ , 6:110–127, January 1933.
  * Ostriker and Peebles [1973] J. P. Ostriker and P. J. E. Peebles.  A numerical study of the stability of flattened galaxies: or, can cold galaxies survive?  _ApJ_ , 186:467–480, December 1973.  doi: 10.1086/152513.
  * Lin and Ishak [2016] Weikang Lin and Mustapha Ishak.  Ultra faint dwarf galaxies: an arena for testing dark matter versus modified gravity.  _JCAP_ , 2016(10):025, October 2016.  doi: 10.1088/1475-7516/2016/10/025.
  * Davis [2000] M. Davis.  The cosmological matter density.  _Phys. Rep._ , 333:147–165, August 2000.  doi: 10.1016/S0370-1573(00)00020-X.
  * Bahcall [1984a] J. N. Bahcall.  K giants and the total amount of matter near the sun.  _ApJ_ , 287:926–944, December 1984a.  doi: 10.1086/162750.
  * Moni Bidin et al. [2012] C. Moni Bidin, G. Carraro, R. A. Méndez, and R. Smith.  Kinematical and Chemical Vertical Structure of the Galactic Thick Disk. II. A Lack of Dark Matter in the Solar Neighborhood.  _ApJ_ , 751(1):30, May 2012.  doi: 10.1088/0004-637X/751/1/30.
  * Garbari et al. [2011] Silvia Garbari, Justin I Read, and George Lake.  Limits on the local dark matter density.  _MNRAS_ , 416(3):2318–2340, 2011.  doi: 10.1051/epjconf/20121901008.
  * Zhang et al. [2013] Lan Zhang, Hans-Walter Rix, Glenn van de Ven, Jo Bovy, Chao Liu, and Gang Zhao.  The gravitational potential near the sun from segue k-dwarf kinematics.  _ApJ_ , 772(2):108, August 2013.  doi: 10.1088/0004-637X/772/2/108.
  * Bienaymé et al. [2014] O. Bienaymé, B. Famaey, A. Siebert, K. C. Freeman, B. K. Gibson, G. Gilmore, E. K. Grebel, J. Bland-Hawthorn, G. Kordopatis, U. Munari, J. F. Navarro, Q. Parker, W. Reid, G. M. Seabroke, A. Siviero, M. Steinmetz, F. Watson, R. F. G. Wyse, and T. Zwitter.  Weighing the local dark matter with RAVE red clump stars.  _A & A_, 571:A92, November 2014.  doi: 10.1051/0004-6361/201424478.
  * Xia et al. [2016] Qiran Xia, Chao Liu, Shude Mao, Yingyi Song, Lan Zhang, RJ Long, Yong Zhang, Yonghui Hou, Yuefei Wang, and Yue Wu.  Determining the local dark matter density with lamost data.  _MNRAS_ , 458(4):3839–3850, 2016.  doi: 10.1093/mnras/stw565.
  * Sivertsson et al. [2018] Sofia Sivertsson, Hamish Silverwood, Justin I Read, Gianfranco Bertone, and Pascal Steger.  The localdark matter density from sdss-segue g-dwarfs.  _MNRAS_ , 478(2):1677–1693, 2018.  doi: 10.1093/mnras/sty977.
  * Buch et al. [2019] Jatan Buch, Shing Chau John Leung, and JiJi Fan.  Using gaia dr2 to constrain local dark matter density and thin dark disk.  _JCAP_ , 2019(04):026, 2019.  doi: 10.1088/1475-7516/2019/04/026.
  * Milgrom [1983] M. Milgrom.  A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis.  _ApJ_ , 270:365–370, July 1983.  doi: 10.1086/161130.
  * Famaey and McGaugh [2012] Benoit Famaey and Stacy McGaugh.  Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions.  _Living rev. relativ._ , 15:10, 2012.  doi: 10.12942/lrr-2012-10.
  * Milgrom [1994] M. Milgrom.  Dynamics with a nonstandard inertia-acceleration relation: An alternative to dark matter in galactic systems.  _Ann. Phys._ , 229(2):384–415, 1994.  ISSN 0003-4916.  doi: https://doi.org/10.1006/aphy.1994.1012.  URL https://www.sciencedirect.com/science/article/pii/S0003491684710128.
  * Lisanti et al. [2019] Mariangela Lisanti, Matthew Moschella, Nadav Joseph Outmezguine, and Oren Slone.  Testing dark matter and modifications to gravity using local milky way observables.  _Phys. Rev. D_ , 100(8):083009, 2019.  doi: 10.1103/PhysRevD.100.083009.
  * Hernandez [2023] X. Hernandez.  Internal kinematics of Gaia DR3 wide binaries: anomalous behaviour in the low acceleration regime.  _MNRAS_ , 525(1):1401–1415, October 2023.  doi: 10.1093/mnras/stad2306.
  * Chae [2023] Kyu-Hyun Chae.  Breakdown of the Newton-Einstein Standard Gravity at Low Acceleration in Internal Dynamics of Wide Binary Stars.  _ApJ_ , 952(2):128, August 2023.  doi: 10.3847/1538-4357/ace101.
  * Pittordis and Sutherland [2023] Charalambos Pittordis and Will Sutherland.  Wide Binaries from GAIA EDR3: preference for GR over MOND?  _OJAp_ , 6:4, February 2023.  doi: 10.21105/astro.2205.02846.
  * Banik et al. [2023] Indranil Banik, Charalambos Pittordis, Will Sutherland, Benoit Famaey, Rodrigo Ibata, Steffen Mieske, and Hongsheng Zhao.  Strong constraints on the gravitational law from Gaia DR3 wide binaries.  _MNRAS_ , November 2023.  doi: 10.1093/mnras/stad3393.
  * Holmberg and Flynn [2004] Johan Holmberg and Chris Flynn.  The local surface density of disc matter mapped by Hipparcos.  _MNRAS_ , 352(2):440–446, August 2004.  doi: 10.1111/j.1365-2966.2004.07931.x.
  * Bovy [2017a] Jo Bovy.  Stellar inventory of the solar neighbourhood using gaia dr1.  _MNRAS_ , 470(2):1360–1387, 2017a.  doi: 10.1093/mnras/stx1277.
  * Marrese, P. M. et al. [2019] Marrese, P. M., Marinoni, S., Fabrizio, M., and Altavilla, G.  Gaia data release 2 - cross-match with external catalogues: algorithms and results.  _A &A_, 621:A144, 2019.  doi: 10.1051/0004-6361/201834142.  URL https://doi.org/10.1051/0004-6361/201834142.
  * Kim and Lépine [2021] Bokyoung Kim and Sebastien Lépine.  Stars in the local galactic thick disc and halo in Gaia EDR3: a catalogue of half a million local main-sequence stars with photometric metallicities.  _MNRAS_ , 510(3):4308–4329, 12 2021.  ISSN 0035-8711.  doi: 10.1093/mnras/stab3671.
  * Pecaut and Mamajek [2013] Mark J Pecaut and Eric E Mamajek.  Intrinsic colors, temperatures, and bolometric corrections of pre-main-sequence stars.  _ApJS_ , 208(1):9, 2013.  doi: 10.1088/0067-0049/208/1/9.
  * Bovy [2017b] Jo Bovy.  Galactic rotation in gaia dr1.  _MNRAS_ , 468(1):L63–L67, 2017b.  doi: 10.1093/mnrasl/slx027.
  * Bekenstein and Milgrom [1984] J. Bekenstein and M. Milgrom.  Does the missing mass problem signal the breakdown of Newtonian gravity?  _ApJ_ , 286:7–14, November 1984.  doi: 10.1086/162570.
  * Navarro et al. [1997] Julio F. Navarro, Carlos S. Frenk, and Simon D. M. White.  A Universal Density Profile from Hierarchical Clustering.  _ApJ_ , 490(2):493–508, December 1997.  doi: 10.1086/304888.
  * Bahcall [1984b] J. N. Bahcall.  The distribution of stars perpendicular to galactic disk.  _ApJ_ , 276:156–168, January 1984b.  doi: 10.1086/161600.
  * Bahcall [1984c] J. N. Bahcall.  Self-consistent determinations of the total amount of matter near the sun.  _ApJ_ , 276:169–181, January 1984c.  doi: 10.1086/161601.
  * Bailer-Jones [2015] Coryn A. L. Bailer-Jones.  Estimating Distances from Parallaxes.  _PASP_ , 127(956):994, October 2015.  doi: 10.1086/683116.
  * Goodman and Weare [2010] Jonathan Goodman and Jonathan Weare.  Ensemble samplers with affine invariance.  _Commun. appl. math. comput. sci._ , 5(1):65–80, 2010.  doi: 10.2140/camcos.2010.5.65.
  * Foreman-Mackey et al. [2013] Daniel Foreman-Mackey, David W Hogg, Dustin Lang, and Jonathan Goodman.  emcee: the mcmc hammer.  _PASP_ , 125(925):306, 2013.  doi: 10.1086/670067.
  * Robert and Wraith [2009] C. P. Robert and D. Wraith.  Computational methods for Bayesian model choice.  In Paul M. Goggans and Chun-Yong Chan, editors, _Bayesian Inference and Maximum Entropy Methods in Science and Engineering: The 29th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering_ , volume 1193 of _American Institute of Physics Conference Series_ , pages 251–262, December 2009.  doi: 10.1063/1.3275622.
  * Newton and Raftery [1994] Michael A. Newton and Adrian E. Raftery.  Approximate bayesian inference with the weighted likelihood bootstrap.  _J. R. Stat. Soc., B (Methodol.)_ , 56(1):3–48, 1994.  ISSN 00359246.  URL http://www.jstor.org/stable/2346025.
  * Malhan et al. [2022] Khyati Malhan, Rodrigo A Ibata, Sanjib Sharma, Benoit Famaey, Michele Bellazzini, Raymond G Carlberg, Richard D’souza, Zhen Yuan, Nicolas F Martin, and Guillaume F Thomas.  The global dynamical atlas of the milky way mergers: Constraints from gaia edr3–based orbits of globular clusters, stellar streams, and satellite galaxies.  _ApJ_ , 926(2):107, 2022.
  * Banik et al. [2017] Nilanjan Banik, Lawrence M. Widrow, and Scott Dodelson.  Galactoseismology and the local density of dark matter.  _MNRAS_ , 464(4):3775–3783, February 2017.  doi: 10.1093/mnras/stw2603.
  * Gentile et al. [2011] Gianfranco Gentile, B Famaey, and WJG de Blok.  Things about mond.  _A &J_, 527:A76, 2011.  doi: 10.1051/0004-6361/201015283.
  * Begeman et al. [1991] K. G. Begeman, A. H. Broeils, and R. H. Sanders.  Extended rotation curves of spiral galaxies : dark haloes and modified dynamics.  _MNRAS_ , 249:523, April 1991.  doi: 10.1093/mnras/249.3.523.
  * Kass and Raftery [1995] Robert E. Kass and Adrian E. Raftery.  Bayes factors.  _J. Am. Stat. Assoc._ , 90(430):773–795, 1995.  doi: 10.1080/01621459.1995.10476572.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572.

##  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:

$\displaystyle p(\mu|x,\sigma)$
$\displaystyle=\frac{p(x|\mu,\sigma)p_{0}(\mu|\sigma)}{E(x,\sigma)}$
$\displaystyle\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

$\displaystyle p(\bm{\mu}|\bm{x},\bm{\sigma})$
$\displaystyle\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

$\displaystyle p(N|\bm{x},\bm{\sigma})$ $\displaystyle=\int
p(N,\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$ $\displaystyle=\int
p(N|\bm{\mu},\bm{x},\bm{\sigma})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$
$\displaystyle=\int
p(N|\bm{\mu})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$
$\displaystyle=\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

$\displaystyle p(\Lambda|\bm{x},\bm{\sigma})$ $\displaystyle=\sum
p(\Lambda,N_{i}|\bm{x},\bm{\sigma})$ $\displaystyle=\sum
p(\Lambda|N_{i},\bm{x},\bm{\sigma})p(N_{i}|\bm{x},\bm{\sigma})$
$\displaystyle=\sum
p(\Lambda|N_{i})\int\mathbb{I}(\bm{\mu})p(\bm{\mu}|\bm{x},\bm{\sigma})\textrm{d}\bm{\mu}$
$\displaystyle=\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:

$\displaystyle
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$:

$\displaystyle\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

$\displaystyle p(\Lambda|\bm{x},\bm{\sigma})$
$\displaystyle\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.