# Exploring Milky Way rotation curves with Gaia DR3: a comparison between
$\Lambda$CDM, MOND, and General Relativistic approaches

Beordo William  Crosta Mariateresa  and Lattanzi Mario Gilberto

###### Abstract

With the release of Gaia DR3, we extend the comparison between dynamical
models for the Milky Way rotation curve initiated in the previous work.
Utilising astrometric and spectro-photometric data for 719143 young disc stars
within $|z|<1$ kpc and up to $R\simeq 19$ kpc, we investigate the accuracy of
MOND and $\Lambda$CDM frameworks in addition to previously studied models,
such as the classical one with a Navarro-Frenk-White dark matter halo and a
general relativistic model. We find that all models, including MOND and
$\Lambda$CDM, are statistically equivalent in representing the observed
rotational velocities. However, $\Lambda$CDM, characterized by an Einasto
density profile and cosmological constraints on its parameters, assigns more
dark matter than the model featuring a Navarro-Frenk-White profile, with the
virial mass estimated at $1.5\text{--}2.5\times 10^{12}\,{\rm M}_{\odot}$ — a
value significantly higher than recent literature estimates. Beyond
$10\text{--}15$ kpc, non-Newtonian/non-baryonic contributions to the rotation
curve are found to become dominant for all models consistently. Our results
suggest the need for further exploration into the role of General Relativity,
dark matter, and alternative theories of gravitational dynamics in shaping
Milky Way’s rotation curve.

##  1 Introduction

With the advent of the new Gaia Data Release 3 (DR3), it is crucial to assess
the accuracy of the dynamical models widely employed in the literature.
Recently, carefully selected rotation curves of the Milky Way (MW) were
presented by [1] for different stellar populations, using data from Gaia DR3.
This work served as a follow-up of [2], which utilised the previous Gaia DR2.
Both studies focused on a direct comparison between a general relativistic
velocity profile derived from the solution of [BG, 3] and a classical
Newtonian model featuring a Navarro-Frenk-White dark matter halo (MWC). In the
present work we extend the analysis of [1] to two state-of-the-art models in
the context of galaxy dynamics: the MOdifyied Newtonian Dynamics [MOND, 4]
paradigm and the $\Lambda$CDM model, where the simulation-motivated Einasto
profile [5, 6] is assumed to describe the CDM density distribution. The
following sections provide a description of the two models adopted and the
data, as well as a discussion and comparison of the results with the MWC and
BG models presented in [1].

##  2 MOND

In the MOND paradigm, the gravitational acceleration is

$\mathbf{g_{\rm MOND}}=\eta\left(\dfrac{g_{\rm N}}{g_{0}}\right)\mathbf{g_{\rm
N}}\,,$ (2.1)

where $\mathbf{g_{\rm N}}$ is the conventional Newtonian acceleration produced
by the baryonic matter alone, while the interpolation function $\eta$ sets the
transition between the Newtonian and the deep MOND regimes through the
acceleration scale $g_{0}$. According to the MOND assumptions, when $g_{\rm
N}\ll g_{0}$ the gravitational acceleration is boosted by a quantity
$\eta\rightarrow\sqrt{g_{0}/g_{\rm N}}$ in order to obtain a flat rotation
curve, while Newtonian dynamics is restored requiring $\eta\rightarrow 1$ when
$g_{\rm N}\gg g_{0}$. Here, we adopt the analytical expression proposed by
[7], namely

$\eta\left(\dfrac{g_{\rm N}}{g_{0}}\right)=\left(1-e^{-\sqrt{g_{\rm
N}/g_{0}}}\right)^{-1}\,.$ (2.2)

This interpolation function has been shown to provide an excellent
representation of the Radial Acceleration Relation (RAR) observed for external
disc galaxies [8, 9].

If we equal the gravitational acceleration $g_{\rm MOND}$ to the centripetal
acceleration $g=V^{2}/R$ of an object orbiting in circular motion with
velocity $V$ at a distance $R$ from the centre of a disc galaxy, the Mondian
rotation velocity results

$V_{\rm MOND}(R,g_{\rm N})=\sqrt{\dfrac{Rg_{\rm N}}{1-e^{-\sqrt{g_{\rm
N}/g_{0}}}}}\,.$ (2.3)

Similarly, the magnitude of the Newtonian acceleration originated by the
distribution of the baryonic matter alone can be written as $g_{\rm N}=V_{\rm
bar}^{2}/R$; therefore the above expression becomes

$V_{\rm MOND}(R,V_{\rm bar})=\dfrac{V_{\rm bar}}{\sqrt{1-e^{-V_{\rm
bar}/\sqrt{Rg_{0}}}}}\,.$ (2.4)

Here, we assume the same density distributions used in [1] for the baryonic
components of the Milky Way (i.e. a Plummer bulge and two Miyamoto-Nagai
discs), so that the MWC and MOND models do not differ where Newtonian gravity
dominates. While in the MWC model of [1] the total rotation curve is given by
adding in quadrature the dark matter halo component (namely $V_{\rm
bar}^{2}+V_{\rm h}^{2}$, where $V_{\rm bar}^{2}=V_{\rm b}^{2}+V_{\rm
td}^{2}+V_{\rm Td}^{2}$), in the MOND model the pure Mondian boost is
represented by the denominator of equation (2.4) and writes explicitly as

$V^{\rm MOND}_{\rm boost}(R,V_{\rm bar})=\sqrt{V_{\rm MOND}^{2}-V_{\rm
bar}^{2}}=V_{\rm bar}\sqrt{\eta(R,V_{\rm bar})-1}\,.$ (2.5)

The free parameters of the baryonic matter distribution share the same prior
distributions of the MWC model [for details see sections 2 and B of 1].

As an additional parameter of the model, we have the acceleration scale
$g_{0}$, which has been constrained to extremely tight values by the observed
RAR of external galaxies [8], namely $g_{0}=(1.20\pm 0.02)\times 10^{-10}$ m
s-2. This scale is supposed to be fixed in the framework of MOND. However,
given the small uncertainty on the parameter, setting a Gaussian prior
$\mathcal{N}(\mu=1.20,\sigma=0.02)\times 10^{-10}$ m s-2 does not affect the
results [10]; therefore, strictly following the Bayesian approach, we prefer
to not fix it and to marginalize over it afterwards.

##  3 $\Lambda$CDM model with Einasto halo profile

We consider the distribution of cold dark matter within the $\Lambda$CDM
scenario to follow the Einasto density profile [5], namely

$\rho_{\rm Einasto}(r)=\rho_{\rm
s}\>\exp\left\\{-\dfrac{2}{\alpha}\left[\left(\dfrac{r}{r_{\rm
s}}\right)^{\alpha}-1\right]\right\\}\;.$ (3.1)

Consistently with the prescriptions of [11], the parameters of the Einasto
profile are written in terms of the halo concentration $C_{200}\equiv
r_{200}/r_{s}$, where the virial radius $r_{200}$ is defined such that the
enclosed average density is 200 times the critical density of the Universe,
i.e. $\rho_{200}=200\rho_{c}=75H_{0}^{2}/(\pi G)$, with $H_{0}=67$ km s-1
Mpc-1 [12, 13]; from here, the rotation velocity and the enclosed halo mass at
the virial radius are then

$V_{200}=10C_{200}r_{s}H_{0}\ ;\
M_{200}=\dfrac{V_{200}^{3}}{10G^{2}H_{0}^{2}}\ .$ (3.2)

With this redefinition, the following boundaries are set: $0<C_{200}<100$,
$10<V_{200}$ [km s-1] $<500$, and $0<\alpha<2$ [11]. Three prior
distributions, coming from $N$-body simulations within the $\Lambda$CDM
cosmology, are then imposed to constrain the parameters:

  * •

Stellar Halo Mass (SHM) relation [14]:

$\dfrac{M_{\star}}{M_{200}}=2N\left[\left(\dfrac{M_{200}}{M_{1}}\right)^{-\beta}+\left(\dfrac{M_{200}}{M_{1}}\right)^{\gamma}\right]^{-1},$
(3.3)

where $M_{\star}=M_{\rm b}+M_{\rm td}+M_{\rm Td}$, $\log(M_{1})=11.59$,
$N=0.0351$, $\beta=1.376$ and $\gamma=0.608$, with scatter $\sigma(\log
M_{\star})=0.15$ dex.

  * •

Halo mass-concentration relation [15]:

$\log(C_{200})=a+b\log(M_{200}/[10^{12}h^{-1}{\rm M}_{\odot}]),$ (3.4)

where $a=0.977$ and $b=-0.130$ assuming the Planck cosmology, with a scatter
of $0.11$ dex.

  * •

Einasto shape parameter as a function of halo mass [13]:

$\alpha=0.0095\lambda^{2}+0.155,$ (3.5)

with a scatter of $0.16$ dex, where
$\log(\lambda)=-0.11+0.146m+0.0138m^{2}+0.00123m^{3}$ and
$m=\log(M_{200}/[10^{12}h^{-1}{\rm M}_{\odot}])$.

Again, the baryonic matter distribution is modelled in the same way as the
MOND and MWC models [see section 2 of 1].

##  4 Data

The data utilized in this study have been carefully selected from Gaia DR3,
following the rigorous criteria outlined in the reference paper [1]. The
resulting sample boasts high-quality astrometric and spectrophotometric
information for a total of 719143 young disc stars located within $|z|<1$ kpc
and spanning from $R=4.5$ to 19 kpc. This includes 241918 O-,B-,A-type stars
(OBA), 475520 Red Giant Branch stars (RGB) with nearly-circular orbits, and
1705 classical Cepheids (DCEP), ensuring a comprehensive representation of
various stellar tracers of the Galactic disc potential. Leveraging this
extensive sample, [1] derived rotation curves of the Milky Way for six
distinct data sets: the 3 ‘pure’ data sets of OBA, DCEP, and RGB stars, the
combined OBA + DCEP and RGB + DCEP samples, and the total sample consisting of
all of the disc stars selected combined (OBA + RGB + DCEP, hereafter ALL).

Here, in order to determine the error bars of the velocity profiles, the
Robust Scatter Estimate (RSE) was adopted as a robust measure of the azimuthal
velocity dispersion of the population in each radial bin, instead of
performing the bootstrapping technique.111The RSE is defined as
$(2\sqrt{2}{\,\rm erf}^{-1}(4/5))^{-1}\approx 0.390152\times(Q90-Q10)$, with
Q90 and Q10 being the $90^{\rm th}$ and $10^{\rm th}$ percentiles of a
distribution, and it coincides with the standard deviation in the case of a
normal distribution. In fact, as discussed in [1], the stellar azimuthal
velocities are considerably dispersed around the median values (the RSE is
typically $>10$ km s-1), therefore smaller error bars would not encompass the
actual variability of the sample. Moreover, with the large amount of stars per
radial bin, a bootstrapped quantity would be much smaller than the individual
uncertainty on the velocity measurements, representing thus a nonphysical
situation.

In the present study, we employ these rotation curves to constrain the models
outlined in sections 2 and 3.

##  5 Results

Table 1 lists, as best-fit estimates, the medians of the posteriors and their
$1\sigma$ credible intervals. Details about the Bayesian analysis can be found
in [1].

MOND model | OBA | DCEP | RGB | OBA $+$ DCEP | RGB $+$ DCEP | ALL
---|---|---|---|---|---|---
$M_{\rm b}$ [$10^{10}$ M⊙] | $1.1^{+0.6}_{-0.6}$ | $1.2^{+0.6}_{-0.6}$ | $1.3^{+0.6}_{-0.6}$ | $1.1^{+0.6}_{-0.6}$ | $1.4^{+0.6}_{-0.6}$ | $1.3^{+0.6}_{-0.6}$
$b_{\rm b}$ [kpc] | $1.0^{+0.9}_{-0.6}$ | $0.8^{+0.8}_{-0.6}$ | $0.8^{+0.8}_{-0.5}$ | $0.9^{+0.9}_{-0.6}$ | $0.8^{+0.7}_{-0.5}$ | $0.8^{+0.8}_{-0.5}$
$M_{\rm td}$ [$10^{10}$ M⊙] | $4.3^{+0.6}_{-0.7}$ | $4.3^{+0.7}_{-0.7}$ | $3.8^{+0.6}_{-0.6}$ | $4.1^{+0.7}_{-0.7}$ | $4.0^{+0.6}_{-0.6}$ | $4.0^{+0.6}_{-0.6}$
$a_{\rm td}$ [kpc] | $5.0^{+0.9}_{-0.8}$ | $5.2^{+0.9}_{-0.8}$ | $4.7^{+0.9}_{-0.9}$ | $4.4^{+0.7}_{-0.6}$ | $4.9^{+0.9}_{-0.8}$ | $4.8^{+0.8}_{-0.8}$
$b_{\rm td}$ [kpc] | $0.3^{+0.4}_{-0.1}$ | $0.4^{+0.5}_{-0.1}$ | $0.4^{+0.5}_{-0.1}$ | $0.3^{+0.3}_{-0.1}$ | $0.4^{+0.5}_{-0.1}$ | $0.3^{+0.5}_{-0.1}$
$M_{\rm Td}$ [$10^{10}$ M⊙] | $4.7^{+0.7}_{-0.7}$ | $4.5^{+0.7}_{-0.8}$ | $4.2^{+0.7}_{-0.8}$ | $4.4^{+0.7}_{-0.7}$ | $4.2^{+0.7}_{-0.7}$ | $4.3^{+0.7}_{-0.8}$
$a_{\rm Td}$ [kpc] | $3.0^{+0.6}_{-0.6}$ | $2.9^{+0.7}_{-0.7}$ | $2.4^{+0.6}_{-0.6}$ | $2.8^{+0.6}_{-0.6}$ | $2.5^{+0.6}_{-0.6}$ | $2.5^{+0.6}_{-0.6}$
$b_{\rm Td}$ [kpc] | $0.9^{+0.8}_{-0.5}$ | $0.8^{+0.9}_{-0.5}$ | $0.4^{+0.7}_{-0.2}$ | $0.7^{+0.7}_{-0.4}$ | $0.5^{+0.7}_{-0.3}$ | $0.5^{+0.7}_{-0.3}$
$g_{0}$ [$10^{-10}$ m s-2] | $1.20^{+0.02}_{-0.02}$ | $1.20^{+0.02}_{-0.02}$ | $1.20^{+0.02}_{-0.02}$ | $1.20^{+0.02}_{-0.02}$ | $1.20^{+0.02}_{-0.02}$ | $1.20^{+0.02}_{-0.02}$
WAIC | $-341\pm 3$ | $-103\pm 3$ | $-345\pm 2$ | $-452\pm 5$ | $-423\pm 5$ | $-426\pm 5$
LOO | $-341\pm 3$ | $-103\pm 3$ | $-345\pm 2$ | $-452\pm 5$ | $-423\pm 5$ | $-426\pm 5$
$\Lambda$CDM model |  |  |  |  |  |
$M_{\rm b}$ [$10^{10}$ M⊙] | $0.8^{+0.6}_{-0.5}$ | $0.9^{+0.6}_{-0.5}$ | $1.1^{+0.6}_{-0.6}$ | $0.8^{+0.6}_{-0.5}$ | $1.1^{+0.6}_{-0.6}$ | $1.0^{+0.6}_{-0.6}$
$b_{\rm b}$ [kpc] | $0.9^{+0.8}_{-0.6}$ | $0.8^{+0.8}_{-0.6}$ | $0.8^{+0.8}_{-0.5}$ | $0.9^{+0.9}_{-0.6}$ | $0.8^{+0.8}_{-0.5}$ | $0.8^{+0.8}_{-0.5}$
$M_{\rm td}$ [$10^{10}$ M⊙] | $3.8^{+0.7}_{-0.7}$ | $3.7^{+0.8}_{-0.8}$ | $3.6^{+0.7}_{-0.7}$ | $4.1^{+0.7}_{-0.7}$ | $3.7^{+0.7}_{-0.7}$ | $3.7^{+0.7}_{-0.7}$
$a_{\rm td}$ [kpc] | $4.8^{+0.9}_{-0.8}$ | $5.1^{+0.9}_{-0.9}$ | $4.7^{+1.0}_{-0.9}$ | $5.0^{+0.8}_{-0.7}$ | $5.0^{+1.0}_{-0.9}$ | $4.9^{+0.9}_{-0.9}$
$b_{\rm td}$ [kpc] | $0.3^{+0.4}_{-0.1}$ | $0.3^{+0.5}_{-0.1}$ | $0.3^{+0.5}_{-0.1}$ | $0.3^{+0.3}_{-0.1}$ | $0.3^{+0.6}_{-0.1}$ | $0.3^{+0.5}_{-0.1}$
$M_{\rm Td}$ [$10^{10}$ M⊙] | $3.9^{+0.9}_{-0.9}$ | $3.7^{+0.9}_{-0.9}$ | $3.8^{+0.9}_{-0.9}$ | $4.2^{+0.8}_{-0.8}$ | $3.7^{+0.8}_{-0.8}$ | $3.8^{+0.8}_{-0.8}$
$a_{\rm Td}$ [kpc] | $2.9^{+0.6}_{-0.7}$ | $2.8^{+0.7}_{-0.7}$ | $2.5^{+0.6}_{-0.6}$ | $3.1^{+0.6}_{-0.7}$ | $2.6^{+0.6}_{-0.6}$ | $2.7^{+0.6}_{-0.6}$
$b_{\rm Td}$ [kpc] | $0.8^{+0.8}_{-0.5}$ | $0.7^{+0.9}_{-0.5}$ | $0.5^{+0.7}_{-0.3}$ | $1.0^{+0.9}_{-0.6}$ | $0.5^{+0.8}_{-0.3}$ | $0.6^{+0.8}_{-0.3}$
$\log(\rho_{\rm s}/[$M⊙ kpc${}^{-3}])$ | $5.8^{+0.3}_{-0.3}$ | $5.9^{+0.3}_{-0.3}$ | $5.8^{+0.3}_{-0.3}$ | $6.0^{+0.3}_{-0.3}$ | $5.9^{+0.3}_{-0.3}$ | $5.9^{+0.3}_{-0.3}$
$r_{\rm s}$ [kpc] | $38^{+20}_{-12}$ | $35^{+16}_{-10}$ | $35^{+17}_{-11}$ | $25^{+10}_{-7}$ | $29^{+11}_{-8}$ | $29^{+11}_{-8}$
$\alpha$ | $0.15^{+0.06}_{-0.04}$ | $0.14^{+0.05}_{-0.04}$ | $0.14^{+0.05}_{-0.04}$ | $0.11^{+0.04}_{-0.03}$ | $0.12^{+0.04}_{-0.03}$ | $0.12^{+0.04}_{-0.03}$
$V_{200}$ [km s-1] | $192^{+27}_{-20}$ | $186^{+22}_{-18}$ | $179^{+20}_{-16}$ | $164^{+12}_{-13}$ | $174^{+13}_{-14}$ | $173^{+13}_{-14}$
$C_{200}$ | $7.5^{+2.4}_{-1.9}$ | $8.0^{+2.5}_{-2.0}$ | $7.7^{+2.5}_{-1.9}$ | $9.8^{+3.0}_{-2.3}$ | $8.8^{+2.5}_{-2.0}$ | $8.8^{+2.6}_{-2.0}$
$\log(M_{200}/[$M${}_{\odot}])$ | $12.39^{+0.17}_{-0.14}$ | $12.35^{+0.15}_{-0.13}$ | $12.30^{+0.14}_{-0.12}$ | $12.19^{+0.09}_{-0.10}$ | $12.26^{+0.10}_{-0.11}$ | $12.26^{+0.09}_{-0.11}$
$r_{200}$ [kpc] | $286^{+40}_{-30}$ | $278^{+33}_{-27}$ | $267^{+30}_{-24}$ | $245^{+17}_{-19}$ | $259^{+20}_{-21}$ | $258^{+19}_{-20}$
WAIC | $-341\pm 4$ | $-103\pm 3$ | $-345\pm 2$ | $-450\pm 5$ | $-425\pm 6$ | $-427\pm 6$
LOO | $-341\pm 4$ | $-103\pm 3$ | $-345\pm 2$ | $-450\pm 5$ | $-425\pm 7$ | $-428\pm 7$

Table 1: Top: estimates of the free parameters of the MOND model, namely, the
medians of the posterior distributions of the Bayesian analysis for each
dataset; the upper and lower bounds (estimated with the $15.87^{\rm th}$ and
$84.13^{\rm th}$ percentiles) enclose their corresponding $1\sigma$ credible
intervals. $M_{\rm b}$, $b_{\rm b}$, $M_{\rm td}$, $a_{\rm td}$, $b_{\rm td}$,
$M_{\rm Td}$, $a_{\rm Td}$, $b_{\rm Td}$, and $g_{0}$ are, respectively, the
mass and the scale length of the bulge, the masses and the scale lengths of
the two discs, and the Mondian acceleration scale. Bottom: estimates of the
free parameters for the $\Lambda$CDM model; besides the parameters of the
baryonic matter distribution, $\rho_{\rm s}$, $r_{s}$, $\alpha$, are the three
parameters defining the Einasto profile, while $V_{200}$, $C_{200}$, $M_{200}$
are respectively the rotation velocity, concentration, and halo mass
calculated at the virial radius $r_{200}$. Log-values of the WAIC and LOO
tests for the Bayesian model comparison are also reported.

The parameters of the baryonic matter distribution are found in agreement
between different datasets and models, even though the $\Lambda$CDM paradigm
tends to assign less mass to the baryonic component, as the values of $M_{\rm
b}$, $M_{\rm td}$, $M_{\rm Td}$ are slightly smaller compared to those
estimated with the MOND and MWC models. Additionally, $\Lambda$CDM
cosmological constraints are observed, as made clear by figure 1:
specifically, the estimated parameters align within 1-$\sigma$ with the
relations for the halo mass-concentration and the Einasto shape parameter, and
within 2-$\sigma$ with the SHM relation.

Figure 1: $\Lambda$CDM constraint relations. Coloured data points are the best
parameter estimates for each stellar sample (see table 1). The shaded gray
regions represent the 1-,2-,3-$\sigma$ ranges around the mean relations.

Rotation curves of the Milky Way for six stellar populations are presented in
figure 2: the results for MOND and $\Lambda$CDM are drawn on top of the
corresponding results of [1] for the BG and MWC models.

Figure 2: The azimuthal velocity profile of the MW as derived from the sample
of disc tracers selected from Gaia DR3. For each dataset, the black starred
symbols represent the median azimuthal velocity at the median distance from
the galactic centre of the stellar population within each of the radial bins
[see tables A.1–A.6 of 1], where the RSE of the velocity distribution defines
the corresponding error bar. The red, blue, green, and yellow curves show the
best-fitting to the BG, MWC, MOND, and $\Lambda$CDM models, respectively. The
filled areas represent the 68 per cent reliability intervals of each rotation
curve; note that for $R\lesssim 4.5$ kpc both the classical and the
relativistic curves are very uncertain because of the lack of data in that
region.

Again, the fits of the four models result statistically equivalent to each
other, having comparable values for the WAIC and LOO tests.

Figure 3 shows the matter density profiles for the four models along with the
only data point considered, which corresponds to the baryonic matter density
observed at the Sun [$\rho_{\rm bar}(R_{\odot})=0.084\pm 0.012$
M${}_{\odot}{\rm pc}^{-3}$, 16].

Figure 3: Density profiles of the MW at $z=0$ for the four models, with their
corresponding 68 per cent confidence intervals; in each panel, the red, blue,
green, and orange solid lines represent the baryonic matter contributions for
the BG, MWC, MOND, and $\Lambda$CDM models respectively. The blue and orange
dashed lines show the total matter for the MWC and $\Lambda$CDM models, while
the dash-dotted lines are the corresponding dark matter contributions. The
vertical grey dashed lines represent the values of $r_{\rm in}$ and $R_{\rm
out}$ of the BG model, while the vertical grey band spans the radial range
covered by the sample. Finally, the black dot represents the local mass
density inferred at the Sun position, i.e. $\rho_{\rm bar}(R_{\odot})=0.084\pm
0.012$ M${}_{\odot}{\rm pc}^{-3}$ from [16].

Estimates of the local baryonic density ($\rho^{\rm\Lambda CDM}_{\rm
bar,\odot}$ and $\rho^{\rm MOND}_{\odot}$, table 2) are consistent with those
in [1], being all around $0.080\,{\rm M}_{\odot}{\rm pc}^{-3}$, as well as
with estimates found in the literature [17, 18, 16]. As for the dark matter
density at the Sun ($\rho^{\rm\Lambda CDM}_{\rm h,\odot}$), we recover again
the recent values reported in the literature [19, 20, 21, 22], being almost
ten times smaller than the local baryonic density: for the full sample,
$\rho^{\rm\Lambda CDM}_{\rm h,\odot}=(0.0090\pm 0.0006)$ M${}_{\odot}{\rm
pc}^{-3}$, which corresponds to $(0.34\pm 0.02)$ GeV cm-3.

Quantity | OBA | DCEP | RGB | OBA $+$ DCEP | RGB $+$ DCEP | ALL
---|---|---|---|---|---|---
$\rho^{\Lambda\rm CDM}_{\rm bar,\odot}$ [${\rm M}_{\odot}{\rm pc}^{-3}$] | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$
$\rho^{\Lambda\rm CDM}_{\rm h,\odot}$ [${\rm M}_{\odot}{\rm pc}^{-3}$] | $0.0099^{+0.0008}_{-0.0008}$ | $0.0098^{+0.0008}_{-0.0008}$ | $0.0088^{+0.0078}_{-0.0007}$ | $0.0085^{+0.0006}_{-0.0006}$ | $0.0091^{+0.0006}_{-0.0006}$ | $0.0090^{+0.0006}_{-0.0006}$
$\rho^{\rm MOND}_{\odot}$ [${\rm M}_{\odot}{\rm pc}^{-3}$] | $0.080^{+0.012}_{-0.012}$ | $0.080^{+0.012}_{-0.012}$ | $0.081^{+0.012}_{-0.012}$ | $0.081^{+0.012}_{-0.012}$ | $0.081^{+0.012}_{-0.012}$ | $0.081^{+0.012}_{-0.012}$
$M^{\Lambda\rm CDM}_{\rm bar}$ [$10^{10}\,{\rm M}_{\odot}$] | $8.49^{+1.10}_{-1.03}$ | $8.31^{+1.11}_{-1.04}$ | $8.51^{+1.06}_{-0.95}$ | $9.23^{+0.95}_{-0.93}$ | $8.43^{+0.96}_{-0.91}$ | $8.49^{+0.96}_{-0.89}$
$M^{\rm MOND}$ [$10^{10}\,{\rm M}_{\odot}$] | $10.12^{+0.33}_{-0.30}$ | $10.05^{+0.48}_{-0.45}$ | $9.32^{+0.27}_{-0.24}$ | $9.59^{+0.19}_{-0.18}$ | $9.59^{+0.21}_{-0.19}$ | $9.56^{+0.21}_{-0.19}$
$M^{\Lambda\rm CDM}_{\rm 200}$ [$10^{12}\,{\rm M}_{\odot}$] | $2.45^{+1.16}_{-0.68}$ | $2.24^{+0.88}_{-0.59}$ | $1.99^{+0.74}_{-0.49}$ | $1.53^{+0.35}_{-0.33}$ | $1.82^{+0.44}_{-0.40}$ | $1.80^{+0.43}_{-0.39}$

Table 2: Estimates of the local baryonic mass density $\rho_{\rm
bar}(R_{\odot})$ for each dataset and both models (respectively
$\rho^{\rm\Lambda CDM}_{\rm bar,\odot}$ and $\rho^{\rm MOND}_{\odot}$) and of
the local dark matter density $\rho^{\rm\Lambda CDM}_{\rm h,\odot}$; the upper
and lower bounds (estimated with the $15.87^{\rm th}$ and $84.13^{\rm th}$
percentiles) enclose their corresponding $1\sigma$ credible intervals.
$M^{\rm\Lambda CDM}_{\rm bar}$ and $M^{\rm MOND}$ are, respectively, the
baryonic mass for the $\Lambda$CDM model and the total mass for the MOND
model; while $M^{\rm\Lambda CDM}_{\rm 200}$ is the estimated virial mass of
the Milky Way in the $\Lambda$CDM scenario. For a direct comparison with the
BG and MWC models, similar estimates can be found in table 2 of [1].

As expected, all density profiles show agreement within the radial range
covered by the data: MWC and $\Lambda$CDM total matter density profiles are
almost coincident, while departing from each other only at very large radii;
so do their dark matter density profiles, but the Einasto profiles of the
$\Lambda$CDM model result larger than the NFW ones both in the inner and outer
parts of the Galaxy; this translates into the dynamical mass being supplied by
more dark matter in the $\Lambda$CDM scenario compared to the case of an NFW
halo without cosmological constraints.

This is confirmed by table 2, where the estimates for the total baryonic and
virial masses are reported (respectively $M^{\Lambda\rm CDM}_{\rm bar}$ and
$M^{\Lambda\rm CDM}_{\rm 200}$). The baryonic mass predicted by the
$\Lambda$CDM scenario is around $8\text{--}9\times 10^{10}\,{\rm M}_{\odot}$,
smaller than the values of $9\text{--}10\times 10^{10}\,{\rm M}_{\odot}$
expected for the MWC model [see table 2 in 1], but still higher compared to
[23], which reports a value of $(5.43\pm 0.57)\times 10^{10}\,{\rm
M}_{\odot}$. The virial mass, ranging in the interval $1.5\text{--}2.5\times
10^{12}\,{\rm M}_{\odot}$, is up to two times the value of $(1.30\pm
0.30)\times 10^{12}\,{\rm M}_{\odot}$ found by [23] and the values in [1] for
the virial mass of the MWC model, significantly higher than those reported by
some authors [24, 19, around $5\text{--}8\times 10^{11}\,{\rm M}_{\odot}$] and
an order of magnitude greater than some very recent values claimed [25, 26,
around $1\text{--}2\times 10^{11}\,{\rm M}_{\odot}$]. These low values from
the literature are the result of a clearly falling rotation curve derived by
the authors beyond 18 kpc; some comments about this surprising finding are
given in section 5.1.

The BG and MOND density profiles are consistent with both the baryonic and
total density profiles of the classical models (MWC and $\Lambda$CDM).
Therefore, no further insight can be gained at this stage without looking at
the amount of mass expected. The total mass predicted by MOND ($M^{\rm MOND}$)
aligns with the total baryonic mass of the MWC model, while being up to 10 per
cent higher than the $\Lambda$CDM estimate. Moreover, this value exceeds the
estimate of $(7.8\pm 0.5)\times 10^{10}\,{\rm M}_{\odot}$ found by [24] in the
MOND framework. Finally, the BG relativistic mass reported in [1] compares
favourably with the value of $1\text{--}2\times 10^{10}\,{\rm M}_{\odot}$
expected for the baryonic mass in $\Lambda$CDM and MOND, calculated in the
region of validity of the BG model.

###  5.1 Comments about the recent claim of a Keplerian rotation curve

Recent claims regarding a Keplerian rotation curve extending up to 30 kpc have
been brought to light [25], a finding seemingly consistent with various recent
studies [27, 24, 26, 28, 19]. Despite the resonance found with these works,
our analysis is restricted by the range of our selected samples, precluding a
direct quantitative comparison. However, several noteworthy distinctions
emerge.

Firstly, within the overlapping range of 10–18 kpc, our rotation curves
exhibit slightly declining profiles, aligning with recent findings that
indicate a pronounced decline only beyond 18–19 kpc (see figure 4).

Figure 4: Zoom-in of the bottom right panel in figure 2, compared with
rotation curves from the recent literature [25, 27, 26, 24, 28].

The profiles of [28] and [26] are found to be the closest match; however, our
rotation curve shows a slight flattening around 15 kpc compared to other
profiles (although statistical limitations of our selected sample start to
become significant at that distance) and exhibits higher velocities in the
inner region close to 5–7 kpc.

Secondly, our approach diverges notably in the methodology employed. While we
refrained from conducting the Jeans analysis, being found negligible in the
DR2 paper [2], we instead implemented an eccentricity selection for the orbits
of RGB stars. This adjustment was aimed at removing the effects of the
asymmetric drift, to match the OBA and DCEP rotation curves. By conducting the
Jeans analysis on our selected sample (as detailed in appendix A), i.e.,
considering the derived circular velocity profile instead of the azimuthal
one, the rotation curve shows a further slight increase within error bars (see
figure 7), as expected. This suggests that the lack of the Jeans analysis in
our procedure is unlikely to be the cause of the discrepancy observed at
around 15 kpc.

Thirdly, our sample selection criteria differed significantly, as we imposed a
stringent requirement of errors on parallaxes smaller than 20%. In contrast,
the referenced literature employed various techniques to extend the measured
rotation curve to 30 kpc, a topic we abstain from delving into here — one of
our rigorous requirements in [1] was to create the most homogeneous sample
ever —, acknowledging the reliability of the methods as asserted by the
authors (see the discussion in the related works).

Furthermore, we delineated our study by confining our analysis to stars within
$|z|<1$ kpc, deviating from the convention of considering a thicker disc
encompassing $|z|<3$ kpc [26]; for instance, a thinner disc selection reflects
into higher velocities, especially between 5–15 kpc from the Galactic Centre
[25, 27].

Additionally, our depiction of error bars on the rotation curves incorporates
the significant dispersion of azimuthal velocity among the stars. As shown in
figure 5, the density plot highlights the azimuthal velocity distribution of
our selected sample.

Figure 5: Same as figure 4 but with a 2D histogram of the full stellar sample
in the $(R,V_{\phi})$-space [see section 3 of 1, for more details about the
data selection].

This stands in contrast to the rotation curves presented by the cited authors,
whose error bars are derived by bootstrapping techniques and typically
encompass only systematic sources of uncertainty, remaining relatively small
compared to the velocity dispersion observed within each radial bin.

Lastly, we adopted quite small radial bins of $0.1$ kpc width (except for the
DCEP sample for which we chose larger bins of $0.5$ kpc) compared to typical
values of 0.5 and 1 kpc used in the literature. As discussed in appendix B,
this of course plays an important role in the derived rotation curve, as the
outer data points are distributed differently (see figure 8). However,
perfectly consistent results of the model parameters are obtained with either
0.5 or 1 kpc radial bins, indicating the robustness of the rotation curve
defined in [1].

These differences in methodology and analysis are potential contributors to
the observed discrepancy in the estimation of the virial mass compared to
recent literature. As a result, our findings suggest that the dynamical mass
of the Galaxy, in the classical context with dark matter, is expected to
remain on the order of $10^{12}\,{\rm M}_{\odot}$.

###  5.2 Contributions to the rotation curve

We can compare the models by making the contributions to the rotation curves
explicit, recalling what was done in section 7 of [1], since different models
and samples contribute to the rotation curve differently. Figure 6, drawn on
top of figure 4 in [1], shows the Newtonian/baryonic counterpart for each
model, alongside the corresponding non-Newtonian/non-baryonic contributions.
The rotation curve of the BG general relativistic model of [1] is split
between an effective Newtonian counterpart and a gravitational dragging
contribution, which is a purely relativistic effect due to the spacetime
geometry. The MWC and $\Lambda$CDM rotation curves are given by the quadratic
summation of the baryonic and halo contributions, while the Mondian boost
responsible for the flat profile is given by equation 2.5. From this plot, it
is clear that the non-Newtonian/non-baryonic contributions onset between
$10\text{--}15$ kpc, becoming dominant beyond 15 kpc. Due to the slightly
higher baryonic mass in MOND, the baryonic component of the rotation curve
tends to be slightly higher compared to the classical models with dark matter.
Similarly, within the framework of $\Lambda$CDM, the dark matter halo
described by the Einasto profile contributes more than the NFW case of the MWC
model, where less dark matter and more baryonic matter are attributed.

Figure 6: Red, blue, green and yellow lines refer to the BG, MWC, MOND, and
$\Lambda$CDM models, respectively. The dashed lines represent the
Newtonian/baryonic counterparts to the rotation curves: $V_{\rm eN}^{\rm BG}$
is the relativistic effective Newtonian velocity for the BG model, $V_{\rm
bar}^{\rm MWC}$ and $V_{\rm bar}^{\rm\Lambda CDM}$ are the velocities
contributed by the baryonic matter for the MWC and $\Lambda$CDM models, and
$V_{\rm bar}^{\rm MOND}$ is the Newtonian velocity contributed by the baryonic
matter for the MOND model. The solid lines show the MWC and $\Lambda$CDM halo
components alone, respectively $V_{\rm h}^{\rm MWC}$ and $V_{\rm
h}^{\rm\Lambda CDM}$, the gravitational dragging contribution for the BG
model, i.e. $V_{\rm drag}^{\rm BG}$, and the Mondian boost $V_{\rm acc}^{\rm
MOND}$. In the bottom right corner of each panel, $z_{\rm eff}$ represents the
effective vertical width of validity for the relativistic disc in the BG
framework.

##  6 Conclusions

In this study, we have undertaken a comprehensive analysis of Galactic
rotation curves utilizing the latest Gaia Data Release 3, in line with seminal
previous results with Gaia DR2 [2]. Considering larger samples from the latest
Gaia Release, we compared the rotation curves within the frameworks of two
prominent dynamical models: MOND and the $\Lambda$CDM paradigm. We extended
the analysis previously presented in [1] with the same set of data, where we
compared again the rotation curves using the BG and MWC models.

Our analysis revealed several key findings. Firstly, we found that both the
MOND and $\Lambda$CDM models provided statistically equivalent fits to the
Gaia DR3 rotation curve data across various stellar populations. This suggests
that all four models are capable of accurately describing the observed
dynamics of the Milky Way, albeit with different underlying assumptions.

The parameters of the baryonic matter distribution are consistent between the
different datasets and models, although the $\Lambda$CDM paradigm tends to
assign slightly less mass to the baryonic component compared to MOND and the
MWC model.

Furthermore, our analysis confirms the previously reported tension between our
derived dynamical mass estimates and those reported in recent literature,
claiming the presence of a Keplerian decline beyond 18-19 kpc. Our estimates
of the virial mass range from $1.5\text{--}2.5\times 10^{12}\,{\rm M}_{\odot}$
in the $\Lambda$CDM framework, which is an order of magnitude higher than
proposed in the literature lately. While our findings align with recent
studies showing near-flat rotation curves within the observed range, we
acknowledge the discrepancy and highlight the need for further investigation
into the methodologies employed.

Additionally, we compared the contributions to the rotation curve from the
baryonic and non-Newtonian components for each model. We found that the non-
Newtonian/non-baryonic contributions become dominant beyond 10–15 kpc, with
MOND predicting a slightly higher baryonic contribution compared to the
classical models with dark matter, and $\Lambda$CDM attributing more dynamical
mass to the dark matter halo described by the Einasto profile compared to the
NFW halo in the MWC model.

Overall, our study underscores the importance of comparing different dynamical
models in understanding the dynamics of galaxies. By leveraging the wealth of
data provided by Gaia DR3, we have gained valuable insights into the structure
and composition of our own Milky Way galaxy, shedding light on the behaviour
and the validity of different theories of gravity. Future work should focus on
resolving the discrepancies observed in the dynamical mass estimates and
further exploring the implications for our understanding of the dark matter
role and galactic structure, especially considering the subtle richness
provided by a general relativistic scenario.

## Data availability

The data sets used in this analysis are presented by [1]. Posterior
distributions or other data underlying this article will be shared on
reasonable request to the corresponding author.

##  Appendix A Jeans analysis of the sample

Under the assumption of axisymmetric and stationary equilibrium, the Jeans
equation in cylindrical coordinates is [29, 30]

$\dfrac{\partial(\nu\langle V_{R}^{2}\rangle)}{\partial
R}+\dfrac{\partial(\nu\langle V_{R}V_{z}\rangle)}{\partial
z}+\nu\left(\dfrac{\langle V_{R}^{2}\rangle-\langle
V_{\phi}^{2}\rangle}{R}+\dfrac{\partial\phi}{\partial R}\right)=0\;,$ (A.1)

where $\nu$ is the radial volume density of the Galaxy, $\phi$ is the
gravitational potential, $V_{R},V_{\phi},V_{z}$ are the radial, azimuthal, and
vertical velocities respectively, and the brackets $\langle\rangle$ represent
quantities averaged over the velocity space. Being the circular rotation
velocity defined as

$V_{c}^{2}=R\dfrac{\partial\phi}{\partial R}\biggr{\rvert}_{z\approx 0}\;,$
(A.2)

the Jeans equation for $z\approx 0$ becomes

$\Biggl{[}\dfrac{\partial(\nu\langle V_{R}^{2}\rangle)}{\partial
R}+\dfrac{\partial(\nu\langle V_{R}V_{z}\rangle)}{\partial
z}+\nu\left(\dfrac{\langle V_{R}^{2}\rangle-\langle
V_{\phi}^{2}\rangle+V_{c}^{2}}{R}\right)\Biggr{]}_{z\approx 0}=0\;.$ (A.3)

Considering the vertical gradient of the cross term $\langle
V_{R}V_{z}\rangle$ negligible, the equation can be written as

$V_{c}^{2}=\left[\langle V_{\phi}^{2}\rangle-\langle
V_{R}^{2}\rangle\left(1+\dfrac{\partial\ln\nu}{\partial\ln
R}+\dfrac{\partial\ln\langle V_{R}^{2}\rangle}{\partial\ln
R}\right)\right]_{z\approx 0}\;,$ (A.4)

where the radial volume density is usually assumed to be $\nu\propto
e^{-R/h_{R}}$. Widely used values for the radial scale length $h_{R}$ range
from 2 to 5 kpc, therefore we adopt a value of 2.5 kpc [31, 25]. The radial
gradient of the averaged squared radial velocity $\langle V_{R}^{2}\rangle$,
i.e. the last term in equation (A.4), can instead be derived directly from the
data set. Fitting $\sqrt{\langle V_{R}^{2}\rangle}$ with an exponential
function, we estimate a scale length of $\approx 26$ kpc, in line with other
studies [26, 28, 19]. The resulting circular velocities for the full sample
are plotted in figure 7. These values typically exceed the azimuthal
velocities by less than 5 per cent and fall well within the error bars, as
expected, given that the orbital eccentricity selection performed in [1] has
already cleared out most of the asymmetric drift.

Figure 7: The top panel shows the circular velocity profile $V_{c}$ derived
from the Jeans equation (A.4) along with the azimuthal velocity $V_{\phi}$ for
the full sample of stars. Here, $V_{c}$ and $V_{\phi}$ share the same error
bars. The bottom panel highlights the relative difference between the two
quantities.

##  Appendix B Dependence on the width of radial bins

The choice of the radial bin size strongly affects the derived rotation curve,
the more the size of the bin is large. This is because some information is
lost and some features of the rotation curve are smoothed out by averaging the
velocities of stars with different true radial coordinates, especially in less
populated bins (i.e. at large radii where data points are crucial in the
dynamical mass determination). For the optimal choice of the bin size, in [1]
we adopted Knuth’s Rule [32], a data-based Bayesian algorithm, yielding to 0.1
kpc bins (0.5 kpc for the DCEP sample). However, in the literature, observed
data are usually grouped in $0.5\text{--}1$ kpc radial bins.

Imposing radial bins of 1 kpc, we find the rotation curve shown in the top
panel of figure 8. The flattening around 15 kpc, although less clear, is still
present; a similar but softer behaviour seems to be found by [26] and [28] at
slightly larger distances and lower velocities.

Figure 8: Same as figure 4 but the azimuthal velocities of the full sample are
averaged over radial bins of size equal to 1 and 0.5 kpc for the left and
right panels respectively.

The bottom panel of the same figure shows the rotation curve obtained with 0.5
kpc radial bins instead, highlighting how data points are differently
distributed with respect to the 1 kpc case. However, if error bars are
appropriately selected, one should obtain consistent results. In fact, when
repeating the fitting process using the rotation curve derived with either 0.5
or 1 kpc bin width for the four models, we obtain estimates of the model
parameters that are perfectly consistent with previous ones. This underlines
the robustness of the rotation curve defined in [1] and used in the present
work. On the contrary, it remains unclear to what extent the results are
influenced by the radial bin size when error bars are derived via
bootstrapping instead, necessitating further verification.

## Acknowledgments

We wish to thank Paola Re Fiorentin and Alessandro Spagna for the selection of
the stellar samples performed in [1]. This work has made use of data products
from: the ESA Gaia mission (gea.esac.esa.int/archive/), funded by national
institutions participating in the Gaia Multilateral Agreement.We are indebted
to the Italian Space Agency (ASI) for their continuing support through
contract 2018-24-HH.0 and its addendum 2018-24-HH.1-2022 to the National
Institute for Astrophysics (INAF).

## References

  * [1] W. Beordo, M. Crosta, M.G. Lattanzi, P. Re Fiorentin and A. Spagna, _Geometry-driven and dark-matter-sustained Milky Way rotation curves with _Gaia_ DR3_, _Monthly Notices of the Royal Astronomical Society_ 529 (2024) 4681.
  * [2] M. Crosta, M. Giammaria, M.G. Lattanzi and E. Poggio, _On testing CDM and geometry-driven Milky Way rotation curve models with Gaia DR2_ , _Monthly Notices of the Royal Astronomical Society_ 496 (2020) 2107.
  * [3] H. Balasin and D. Grumiller, _NON-NEWTONIAN BEHAVIOR IN WEAK FIELD GENERAL RELATIVITY FOR EXTENDED ROTATING SOURCES_ , _Int. J. Mod. Phys. D_ 17 (2008) 475.
  * [4] M. Milgrom, _A modification of the Newtonian dynamics - Implications for galaxies._ , _The Astrophysical Journal_ 270 (1983) 371.
  * [5] J. Einasto, _On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters_ , _Trudy Astrofizicheskogo Instituta Alma-Ata_ 5 (1965) 87.
  * [6] J.F. Navarro, E. Hayashi, C. Power, A.R. Jenkins, C.S. Frenk, S.D.M. White et al., _The inner structure of $\Lambda$CDM haloes - III. Universality and asymptotic slopes_, _MNRAS_ 349 (2004) 1039 [astro-ph/0311231].
  * [7] S.S. McGaugh, _Milky Way Mass Models and MOND_ , _ApJ_ 683 (2008) 137.
  * [8] F. Lelli, S.S. McGaugh, J.M. Schombert and M.S. Pawlowski, _One Law to Rule Them All: The Radial Acceleration Relation of Galaxies_ , _ApJ_ 836 (2017) 152.
  * [9] S.S. McGaugh, F. Lelli and J.M. Schombert, _Radial Acceleration Relation in Rotationally Supported Galaxies_ , _Phys. Rev. Lett._ 117 (2016) 201101.
  * [10] P. Li, F. Lelli, S. McGaugh and J. Schombert, _Fitting the radial acceleration relation to individual SPARC galaxies_ , _A &A_ 615 (2018) A3.
  * [11] P. Li, F. Lelli, S.S. McGaugh, N. Starkman and J.M. Schombert, _A constant characteristic volume density of dark matter haloes from SPARC rotation curve fits_ , _Monthly Notices of the Royal Astronomical Society_ 482 (2019) 5106.
  * [12] Planck Collaboration, P.A.R. Ade, N. Aghanim, C. Armitage-Caplan, M. Arnaud, M. Ashdown et al., __Planck_ 2013 results. XVI. Cosmological parameters_, _A &A_ 571 (2014) A16.
  * [13] A.A. Dutton and A.V. Macciò, _Cold dark matter haloes in the Planck era: Evolution of structural parameters for Einasto and NFW profiles_ , _Monthly Notices of the Royal Astronomical Society_ 441 (2014) 3359.
  * [14] B.P. Moster, T. Naab and S.D.M. White, _Galactic star formation and accretion histories from matching galaxies to dark matter haloes_ , _Monthly Notices of the Royal Astronomical Society_ 428 (2013) 3121.
  * [15] A.V. Macciò, A.A. Dutton and F.C. Van Den Bosch, _Concentration, spin and shape of dark matter haloes as a function of the cosmological model:_WMAP_ 1, _WMAP_ 3 and _WMAP_ 5 results_, _Monthly Notices of the Royal Astronomical Society_ 391 (2008) 1940.
  * [16] C.F. McKee, A. Parravano and D.J. Hollenbach, _STARS, GAS, AND DARK MATTER IN THE SOLAR NEIGHBORHOOD_ , _ApJ_ 814 (2015) 13.
  * [17] S. Garbari, C. Liu, J.I. Read and G. Lake, _A new determination of the local dark matter density from the kinematics of K dwarfs_ , _MNRAS_ 425 (2012) 1445 [1206.0015].
  * [18] O. Bienaymé, B. Famaey, A. Siebert, K.C. Freeman, B.K. Gibson, G. Gilmore et al., _Weighing the local dark matter with RAVE red clump stars_ , _A &A_ 571 (2014) A92 [1406.6896].
  * [19] A.-C. Eilers, D.W. Hogg, H.-W. Rix and M.K. Ness, _The Circular Velocity Curve of the Milky Way from 5 to 25 kpc_ , _The Astrophysical Journal_ 871 (2019) 120.
  * [20] M. Cautun, A. Benítez-Llambay, A.J. Deason, C.S. Frenk, A. Fattahi, F.A. Gómez et al., _The milky way total mass profile as inferred from Gaia DR2_ , _MNRAS_ 494 (2020) 4291 [1911.04557].
  * [21] A. Widmark, C.F.P. Laporte, P.F. de Salas and G. Monari, _Weighing the Galactic disk using phase-space spirals: II. Most stringent constraints on a thin dark disk using _Gaia_ EDR3_, _A &A_ 653 (2021) A86.
  * [22] J. Wang, F. Hammer and Y. Yang, _Milky Way total mass derived by rotation curve and globular cluster kinematics from Gaia EDR3_ , _Monthly Notices of the Royal Astronomical Society_ 510 (2022) 2242.
  * [23] P.J. McMillan, _The mass distribution and gravitational potential of the Milky Way_ , _Monthly Notices of the Royal Astronomical Society_ 465 (2017) 76.
  * [24] F.S. Labini, Ž. Chrobáková, R. Capuzzo-Dolcetta and M. López-Corredoira, _Mass Models of the Milky Way and Estimation of Its Mass from the Gaia DR3 Data Set_ , _ApJ_ 945 (2023) 3.
  * [25] Y. Jiao, F. Hammer, H. Wang, J. Wang, P. Amram, L. Chemin et al., _Detection of the Keplerian decline in the Milky Way rotation curve_ , _A &A_ 678 (2023) A208.
  * [26] X. Ou, A.-C. Eilers, L. Necib and A. Frebel, _The dark matter profile of the Milky Way inferred from its circular velocity curve_ , _Monthly Notices of the Royal Astronomical Society_ 528 (2024) 693.
  * [27] H.-F. Wang, Ž. Chrobáková, M. López-Corredoira and F. Sylos Labini, _Mapping the Milky Way Disk with Gaia DR3: 3D Extended Kinematic Maps and Rotation Curve to $\approx$30 kpc_, _ApJ_ 942 (2023) 12.
  * [28] Y. Zhou, X. Li, Y. Huang and H. Zhang, _The Circular Velocity Curve of the Milky Way from 5–25 kpc Using Luminous Red Giant Branch Stars_ , _ApJ_ 946 (2023) 73.
  * [29] J.H. Jeans, _On the theory of star-streaming and the structure of the universe_ , _MNRAS_ 76 (1915) 70.
  * [30] J. Binney and S. Tremaine, _Galactic Dynamics_ , Princeton Series in Astrophysics, Princeton University Press, Princeton, 2nd ed. (2008).
  * [31] M. Jurić, Ž. Ivezić, A. Brooks, R.H. Lupton, D. Schlegel, D. Finkbeiner et al., _The Milky Way Tomography with SDSS. I. Stellar Number Density Distribution_ , _ApJ_ 673 (2008) 864 [astro-ph/0510520].
  * [32] K.H. Knuth, _Optimal Data-Based Binning for Histograms_ , _arXiv e-prints_ (2006) physics/0605197.