11institutetext:  NICPB, Rävala 10, Tallinn 10143, Estonia 11email:
<EMAIL_ADDRESS>22institutetext: Tallinn University of Technology,
Ehitajate tee 5, Tallinn 19086, Estonia 33institutetext: Tartu Observatory,
University of Tartu, Observatooriumi 1, Tõravere 61602, Estonia 33email:
<EMAIL_ADDRESS>44institutetext: Instituto de Astrofísica de
Canarias, C/Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain 55institutetext:
Universidad de La Laguna, Dpto. Astrofísica, Avenida Astrofísico Francisco
Sánchez, 38206 La Laguna, Tenerife, Spain

# A Bayesian estimation of the Milky Way’s circular velocity curve using Gaia
DR3

Sven Põder 1122 María Benito 3311 Joosep Pata 11 Rain Kipper 33 Heleri Ramler
33 Gert Hütsi 11 Indrek Kolka 33 Guillaume F. Thomas 4455

###### Abstract

Aims. Our goal is to calculate the circular velocity curve of the Milky Way,
along with corresponding uncertainties that quantify various sources of
systematic uncertainty in a self-consistent manner.

Methods. The observed rotational velocities are described as circular
velocities minus the asymmetric drift. The latter is described by the radial
axisymmetric Jeans equation. We thus reconstruct the circular velocity curve
between Galactocentric distances from 5 kpc to 14 kpc using a Bayesian
inference approach. The estimated error bars quantify uncertainties in the
Sun’s Galactocentric distance and the spatial-kinematic morphology of the
tracer stars. As tracers, we used a sample of roughly 0.6 million stars on the
red giant branch stars with six-dimensional phase-space coordinates from Gaia
data release 3 (DR3). More than 99% of the sample is confined to a quarter of
the stellar disc with mean radial, rotational, and vertical velocity
dispersions of $(35\pm 18)\,\rm km/s$, $(25\pm 13)\,\rm km/s$, and $(19\pm
9)\,\rm km/s$, respectively.

Results. We find a circular velocity curve with a slope of $0.4\pm 0.6\,\rm
km/s/kpc$, which is consistent with a flat curve within the uncertainties. We
further estimate a circular velocity at the Sun’s position of
$v_{c}(R_{0})=233\pm 7\,\rm km/s$ and that a region in the Sun’s vicinity,
characterised by a physical length scale of $\sim 1\,\rm kpc$, moves with a
bulk motion of $V_{LSR}=7\pm 7\,\rm km/s$. Finally, we estimate that the dark
matter (DM) mass within 14 kpc is $\log_{10}M_{\rm DM}(R<14\,{\rm kpc})/{\rm
M_{\odot}}=\left(11.2^{+2.0}_{-2.3}\right)$ and the local spherically averaged
DM density is $\rho_{\rm DM}(R_{0})=\left(0.41^{+0.10}_{-0.09}\right)\,{\rm
GeV/cm^{3}}=\left(0.011^{+0.003}_{-0.002}\right)\,{\rm M_{\odot}/pc^{3}}$. In
addition, the effect of biased distance estimates on our results is assessed.

###### Key Words.:

Galaxy: kinematics and dynamics - Galaxy: disk - Methods: statistical

##  1 Introduction

The rotation of stars and gas in galactic discs has been extensively used as a
kinematical tracer of matter distribution of external galaxies and our own
Galaxy, the Milky Way (MW) (Kuzmin 1970). Several recent studies (Mróz et al.
2019; Eilers et al. 2019; Chrobáková et al. 2020; Ablimit et al. 2020; Khanna
et al. 2023; Gaia Collaboration et al. 2023; Wang et al. 2022; Zhou et al.
2023) have measured the stellar disc rotation in the MW using stellar data
from the Gaia satellite (Gaia Collaboration et al. 2016). These studies differ
in the samples of stars used as a tracer and/or in the methodology and
assumptions employed. Moreover, some of the cited studies provided rotational
(azimuthal) velocities, whereas the others presented circular ones. The former
are direct measurements and no underlying assumptions are made with respect to
the shape or time dependence of the MW’s gravitational potential. On the other
hand, circular velocities assume a stationary gravitational potential that
exhibits axial symmetry. Moreover, these are the velocities that should be
used to derive the total or dynamical matter distribution. The modelling
assumptions can therefore bias the determination of the total and dark matter
content in our Galaxy.

The amount of phase space data currently available is large, and statistics is
generally not the limiting factor for studies of the dynamics of the MW
stellar disc. The limiting factor is instead systematic errors, such as the
Sun’s Galactocentric distance or the adoption of incorrect modelling
assumptions. In this paper, we present a Bayesian inference approach to
estimate the circular velocity curve of our Galaxy that allows the
straightforward incorporation of systematic and statistical sources of
uncertainty, which are both treated as nuisance parameters. In this way, we
provide, for the first time, a circular velocity curve with errorbars that
self-consistently include uncertainties in the Sun’s Galactocentric distance
and in the spatial-kinematic structure of the stellar disc. Therefore, our
uncertain knowledge about astrophysical parameters is propagated through
Bayes’ theorem to the determination of the circular velocity curve and,
subsequently, to the estimation of the dark matter density profile in our
Galaxy. Taking into account the uncertainties about how dark matter is
distributed in the MW is essential for interpreting the results of dark matter
particle searches (see e.g. Benito et al. 2017).

The Bayesian inference approach presented here is a flexible method that
models the observed rotational or azimuthal velocity at a given Galactocentric
distance as the difference between the circular velocity and the asymmetric
drift component. The latter velocity component is obtained from the
stationary, axisymmetric radial Jeans equation under the assumption of
symmetry above and below the Galactic plane. Observed and modelled azimuthal
velocities are then compared by means of the Bayes theorem. The paper is
divided as follows. Section 2 describes the data; Section 3 presents the
Bayesian methodology. The results are presented in 4, and we conclude in
Section 5.

##  2 Data

###  2.1 Red giant branch stars

We used astrometric data and radial velocities from Gaia data release 3 (DR3)
for $665\,660$ stars in the red giant branch (RGB). These stars are old and
have relatively large velocity dispersion. Thus, they are less susceptible to
perturbations. They are also bright enough to have measured radial velocities
out to large distances. Specifically, we used the same sample of almost six
million RGB stars as in Gaia Collaboration et al. (2023) to which we performed
additional spatial and kinematic cuts.

In the following we briefly describe how the RGB sample was obtained and we
refer readers to the original paper Gaia Collaboration et al. (2023) for a
thorough description. Red giants are identified as stars with effective
temperatures between $3000\,{\rm K}$ and $5500\,\rm K$ and surface gravity
satisfying the condition $\log g<3$. Both stellar parameters are provided as
data products in Gaia DR3 (Andrae, R. et al., 2023). Using these first
selection criteria, we obtained 11 576 957 sources. We then selected RGB stars
with good astrometric data as measured by the fidelity parameter $f_{a}$ given
in Rybizki et al. (2022) and we removed stars with $f_{a}\leq 0.5$, thus
remaining with 6 586 329 stars. After this, we performed a cut in height above
and below the Galactic plane, $|z|<0.2\,\rm kpc$ (which removes ca. 4.5M
stars), a cut in Galactocentric distances
$$5\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$<R<$14\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$
(removing ca. 84k stars), and a cut in heliocentric velocity
$|\textbf{v}-\textbf{v}_{\odot}|<210\,\rm km/s$.111$\textbf{v}_{\odot}$ is
defined in Eq. (1). The latter cut removed roughly 28k stars. The $z$ and
velocity cuts are applied to remove halo stars (Helmi et al. 2018; Thomas et
al. 2019). The cut in height also avoids large density and velocity gradients
in the z-coordinate, thus making the derivatives with respect to $z$ in the
Jeans equation negligible (see Eq. (3)). We note that the scale height of the
thin disc is roughly 250 pc (Bland-Hawthorn & Gerhard 2016).

Figure 1: Bias in different distance estimates. Left: GSP-Phot heliocentric
distances as a function of parallax inverse. The colour bar shows the mean
parallax quality inside each pixel. The quality cut used in this paper
corresponds to $\log_{10}\left(\varpi/\sigma_{\varpi}\right)>1.301$. Right:
same but for photogeo heliocentric distances by Bailer-Jones et al. (2021).

We used GSP-Phot distances (Andrae, R. et al., 2023) instead of the inverse of
the parallax due to noisy parallax measurements. As shown in the left panel of
Fig. 1, the GSP-Phot distances are significantly underestimated at large
distances from the Sun (see also Andrae, R. et al. 2023; Fouesneau, M. et al.
2023). This would lead to artificially lower circular velocities and thus to a
steeper slope of the estimated circular velocity curve. To reduce this bias,
we imposed a tight constraint on the quality of the parallax measurements,
namely $\varpi/\sigma_{\varpi}>20$. This cut-off removed approximately 1.3M
stars and suppressed the systematic underestimation of the distances, in fact
leading to a slight overestimation. To assess the dependence of our results on
inaccurate distances, we further performed the same analysis using a less
stringent quality cut on parallaxes and ’photogeo’ distances from Bailer-Jones
et al. (2021) (henceforth referred to as BJ21). Photogeo distances suffer from
underestimation when including measurements with
$\varpi/\sigma_{\varpi}\lesssim 20$ (see the right panel of Fig. 1). In the
end, we are left with a final sample of $665\,660$ stars which is shown in
Figure 2.

Figure 2: Spatial distribution of the final red giants sample. Left:
distribution in the Galactic longitude ($\ell$) and latitude ($b$) plane. The
colour bar shows the number density of stars per pixel. Each pixel has a size
of 1 degree in both latitude and longitude. Right: same but projected into the
Galactic plane. Each pixel has a size of 0.05 kpc and 0.06 kpc in the x and y
coordinates, respectively. In this figure, the Galactic centre is located at
(0, 0), the Sun is located at $(-8.277,0)$ and the rotation of the Galaxy is
clockwise. We added dashed, black circles at 5 kpc,
$R_{0}=$8.277\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$, 12 kpc and 14 kpc to
ease visualisation. For the transformation to Galactocentric coordinades, we
adopt $R_{0}=8.277\,\rm kpc$, $z_{0}=25\,\rm pc$ and ${\bf
v}_{\odot}=(11.1,251.5,8.59)\,\rm km/s$.

Figure 3: Distribution of the mean observed rotational velocity $v_{\phi}$
inside x-y pixels for our final sample of RGB stars (left). The size of the
bin is 150 pc in both x and y coordinates, respectively. We added dashed,
black circles at 5 kpc,
$R_{0}=$8.277\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$$, 12 and 14 kpc to ease
visualisation. Square root of the radial and azimuthal diagonal components of
the velocity dispersion tensor in radial bins of 1 kpc (right).

###  2.2 Galactocentric frame

####  2.2.1 Galactocentric transformation

We transformed RGB stars from a heliocentric to a Galactocentric reference
frame. We used a right-handed Galactocentric coordinate system $(x,y,z)$, with
the Sun located at negative $x$, $y$ pointing in the direction of the Galactic
rotation and $z$ towards the North Galactic Pole. In order to define this new
frame and perform the transformation correctly, we assumed solar orbital
parameters ($R_{0}$, $z_{0}$, $U_{\odot}$, $V_{\odot}$, $W_{\odot}$) from
contemporary literature. First, we treated the Galactocentric distance $R_{0}$
as a nuisance parameter of the analysis in order to account for uncertainties
in its determination. In particular, we used a uniform prior range
$R_{0}\in[7.8-8.5]\,\rm kpc$ that encompasses recent estimates with their
corresponding uncertainties (Do et al. 2019; GRAVITY Collaboration et al.
2021; Abuter et al. 2018; GRAVITY Collaboration et al. 2019, 2020, 2022; Leung
et al. 2022). Our intention was not to constrain the value of $R_{0}$, but
rather to show how the uncertainty in this parameter, which is encoded in the
prior, propagates into the circular velocity curve. For this reason, we
remained agnostic about the actual value of $R_{0}$ and adopted a uniform
distribution that encompasses the most recent $R_{0}$ estimations within
2$\sigma$ uncertainty. Second, for the height of the Sun over the Galactic
plane $z_{0}$, we assumed a value of 25 pc (Jurić et al. 2008)222By adopting
an alternative $Z_{0}$ value of 0 pc, the change in the central circular
velocities is smaller than 1%.. The transformation from spherical coordinates
in ICRS to a Galactocentric Cartesian setting was done as described in Johnson
& Soderblom (1987) and Hobbs et al. (2018).

Finally, the radial velocity measurements from Gaia also allow us to construct
the full 3D space velocities and then transform them to a new frame by using
the Sun’s Galactocentric velocity. Under the assumption that Sg A∗ is at rest
at the Galactic centre, the y and z components of the total solar velocity
vector are derived from the proper motion of Sg A∗, as measured in Reid &
Brunthaler (2020), in combination with the adopted value of $R_{0}$. On the
other hand, for the x-component of the velocity, we adopted the value from
Schönrich et al. (2010). We have not corrected this value for the offset in
radial velocity between the radio-to-infrared reference frames determined by
the GRAVITY collaboration (GRAVITY Collaboration et al. 2019, 2020, 2021,
2022) as suggested in Drimmel & Poggio (2018). There are many possible sources
for this systematic offset and, in any case, it is compatible with zero at the
2$\sigma$ level (GRAVITY Collaboration et al. 2022). In this way, we obtain
the following vector

${\bf v}_{\odot}=\begin{bmatrix}U_{\odot}\\\ V_{\odot}\\\
W_{\odot}\end{bmatrix}=\begin{bmatrix}11.1\\\
251.5\times\left(\frac{R_{0}}{$8.277\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}\right)\\\
8.59\times\left(\frac{R_{0}}{$8.277\text{\,}\mathrm{k}\mathrm{p}\mathrm{c}$}\right)\end{bmatrix}\rm
km/s,$ (1)

where $U_{\odot}$, $V_{\odot}$, $W_{\odot}$ correspond to the velocity
components in the Galactocentric x, y, and z-directions, respectively. As
$U_{\odot}$ corresponds to the radial motion of the Sun towards the Galactic
centre, we implicitly assumed that the LSR has no such motion.

Using Gaia measurements for right ascension, declination, and radial velocity,
and the GSP-Phot distances with a quality parallax cut of
$\varpi/\sigma_{\varpi}>20$, we transformed the proper motions and radial
velocities first to Cartesian velocities in a similar way as was done in
Johnson & Soderblom (1987) and Hobbs et al. (2018). Finally, we switched to
Galactocentric cylindrical coordinates $(R,\phi,z,v_{R},v_{\phi},v_{z})$.

The left panel of figure 3 shows the mean observed rotational velocities in
the Galactic plane. The right panel of the same figure depicts the mean
axisymmetric radial and azimuthal velocity dispersions. Figure 4 shows the
mean azimuthal and radial velocities. In both figures uncertainties are
calculated by bootstrap resampling and are given by half the interval between
the 16th and 84th percentiles of the corresponding distribution. As shown in
the right panel of the last figure, the bulk motion of the stars in the radial
direction exhibits an oscillatory pattern with an amplitude of roughly 5 km/s.
This was already reported in Gaia DR2 (Gaia Collaboration et al., 2018) and
might be a kinematic signature of the spiral arms or the result of interaction
with a perturber. We leave for future work a careful study of the origin of
this intriguing oscillation in the radial velocities.

Figure 4: Mean azimuthal (left) and radial (right) velocities as a function of
Galactocentric distance $R$. As in figure 3, errorbars are calculated via
bootstrap (see text for details).

####  2.2.2 Binning

We treated the Sun’s Galactocentric distance as a free parameter in our
analysis. Varying $R_{0}$ translates into a variation of the $R$ coordinate of
the RGB stars. For this reason, we distributed the RGB sample in bins defined
in the adimensional coordinate $x=R/R_{0}$. This mitigates the shift of the
red giants’ $R$ coordinate when varying $R_{0}$. Thus, for different $R_{0}$
values, a given $x$ bin contains approximately the same stars (Benito et al.
2019). Furthermore, we note that by marginalising over the azimuthal
coordinate $\phi$, the rotational velocity in the Galactic disc is treated as
a purely radially dependent observable.

In total, we defined eight bins from $x=5/8.277$ to $x=14/8.277$ with a step
of $\Delta x=1/8.277$. Observed rotational (azimuthal) velocities inside each
bin approximately follow a Gaussian distribution as shown in figure 5. In this
figure, we plot the distribution of observed velocities inside two bins: the
bin where the distribution deviates most from a Gaussian and a randomly
selected bin. The rotational velocity inside each bin is defined by the mean.
The associated uncertainties are calculated by bootstrap resamplings and are
given by half the interval between the 16th and 84th percentiles of the
velocity distributions.

Figure 5: Distribution of $v_{\phi}$ inside the last radial bin where the
distribution of values deviates most from a Gaussian distribution (orange) and
a randomly selected bin (blue). The dashed lines depict the best-fitted
Gaussian of the rotational velocity inside each particular bin.

##  3 Methodology

###  3.1 Axisymmetric kinematic model

Inside each radial bin, we modeled the mean rotational or azimuthal velocity
as

$v_{\phi{{\rm model}}}=v_{c}-v_{a},$ (2)

where $v_{c}$ is the circular velocity or the velocity of a star moving in a
circular orbit and $v_{a}$ is the asymmetric drift. The latter accounts for
the diffusion of stars in phase-space as the stars orbit the Galaxy and
streaming or bulk motions inside the disc. Neglecting large scale non-
axisymmetric features, the asymmetric drift component can be obtained from the
radial Jeans equation under the assumption that the MW is in a steady-state
and has an axisymmetric gravitational potential (Binney & Tremaine 2008). If
we further expect the density distribution to be symmetric with respect to the
Galactic plane, the Jeans equation takes the following form

$\frac{R}{\nu}\frac{\partial\left(\nu\overline{v_{R}^{2}}\right)}{\partial
R}+R\frac{\partial\overline{v_{R}v_{z}}}{\partial
z}+\overline{v_{R}^{2}}-\overline{v_{\phi}^{2}}+v_{c}^{2}=0.$ (3)

Substituting
$\overline{v_{\phi}^{2}}=\sigma_{\phi}^{*2}+(\overline{v_{\phi}})^{2}$
iiiiiiiiiTo avoid confusion between uncertainties and the components of the
velocity-dispersion tensor, we add the superscript $*$ to the latter as in
Gaia Collaboration et al. (2023). and assuming
$\overline{v_{R}^{2}}=\sigma_{R}^{*2}$, we further obtain

$\sigma_{\phi}^{*2}-\sigma_{R}^{*2}-\frac{R}{\nu}\frac{\partial\left(\nu\sigma_{R}^{*2}\right)}{\partial
R}-R\frac{\partial\overline{v_{R}v_{z}}}{\partial
z}=v_{c}^{2}-(\overline{v_{\phi}})^{2}.$ (4)

We have checked the latter assumption explicitly and observed that the
measured values of $\overline{v_{R}^{2}}$ are at least two orders of magnitude
larger than the measured $\overline{v_{R}}^{2}$.

In traditional approaches, circular velocities in radial bins are calculated
by plugging numbers into Eq. (4). Circular velocity values between bins are
correlated as the radial dependence of the density profile and radial velocity
dispersion are described by exponential functions. In the proposed new
approach, we add Eq. (2) to transform the problem into an inference procedure.
In this way, we can introduce free parameters such as scale lengths of the
exponential functions and/or the Sun’s Galactocentric distance $R_{0}$. The
fitting procedure introduces an additional correlation between the bins.
Nonetheless, we found that if we fix all free parameters of the analysis (i.e.
scale lengths and $R_{0}$), the traditional and new approaches give the same
circular velocity curve.

By expanding the difference of squares on the right hand side and introducing
equation (2), we arrive to the following expression for the asymmetric drift

$v_{a}=\frac{\sigma_{R}^{*2}}{v_{c}+\overline{v_{\phi}}}\left[\frac{\sigma_{\phi}^{*2}}{\sigma_{R}^{*2}}-1-\frac{\partial\ln\nu}{\partial\ln
R}-\frac{\partial\ln(\sigma_{R}^{*2})}{\partial\ln
R}-\frac{R}{\sigma_{R}^{*2}}\frac{\partial\overline{v_{R}v_{z}}}{\partial
z}\right].$ (5)

The sum $v_{c}+\overline{v_{\phi}}$ in the denominator is often approximated
as $\approx 2v_{c}$ (Binney & Tremaine 2008). Nonetheless, we leave it as it
is and estimate $\overline{v_{\phi}}$ as the mean rotational velocity inside
each bin. In addition to this, the diagonal components of the velocity-
dispersion tensor $\sigma_{\phi}^{*2}$ and $\sigma_{R}^{*2}$ are also
calculated directly from the data, and they correspond to the variance of the
azimuthal and radial velocity in the bin respectively. For the 3rd component
inside the brackets of equation (5), the number density distribution $\nu$ is
described by an exponential profile, namely $\nu\propto\exp{(-R/h_{R})}$ with
$h_{R}$ the disc scale radius. Notice that in the 4th component we describe
$\sigma_{R}^{*}$ as $\sigma_{R}^{*}\propto\exp(-R/h_{\sigma})$, where
$h_{\sigma}$ is the scale length of the radial velocity dispersion. Finally,
after taking the derivatives with respect to $\ln R$ in (5), we are left with
the following equation

$v_{a}=\frac{\sigma_{R}^{*2}}{v_{c}+\overline{v_{\phi}}}\left[\frac{\sigma_{\phi}^{*2}}{\sigma_{R}^{*2}}-1+R\left(\frac{1}{h_{r}}+\frac{2}{h_{\sigma}}\right)\right].$
(6)

We neglect the last term of (5) in our axisymmetric treatment of the rotation
curve derivation as $\overline{v_{R}v_{z}}\approx 0$. This is motivated
because the radial and vertical motions are expected to decouple for circular
orbits near the disc when the velocity ellipsoid is aligned with the Galactic
plane (Bovy 2023). In reality, however, this is not necessarily true (e.g.
some general models are provided by Tempel & Tenjes 2006; Kipper et al. 2016).
In any case, our cut in z-coordinate $|z|<0.2\,\rm kpc$ minimises gradients in
the vertical direction. Furthermore, the inclusion of this term changes the
final circular velocity at the percent level, as shown in Eilers et al.
(2019).

###  3.2 Circular velocity fitting

We used the axisymmetric kinematic model described in the previous section to
derive the circular velocity $v_{c,j}$ in each radial bin $j$. We approached
this as a Bayesian inference problem, wherein we used a Markov chain Monte
Carlo (MCMC) algorithm to sample the posterior probability of our model
parameters $\theta$, namely the circular velocities $\\{v_{c,j}\\}$, $h_{R}$,
$h_{\Theta}$ and $R_{0}$. According to Bayes’ theorem, the posterior
distribution of a set of model parameters $\theta$ given a particular set of
data $D$ can be defined as

$p(\theta|{D})=\frac{p({D}|\theta)p(\theta)}{p(D)},$ (7)

where $p(\theta)$ is the prior distribution function that contains a priori
knowledge about the parameters, $p(D)$ is the Bayesian evidence which is an
irrelevant normalisation constant in this context. Moreover, the likelihood
function takes the form

$p(D|\theta)=-\prod_{j}\left[\frac{1}{\sqrt{2\pi\sigma_{\overline{v_{\phi,j}}}^{2}}}\exp\left(\frac{\left(\overline{v_{\phi,j}}-v_{\phi{{\rm
model},j}}(\theta)\right)^{2}}{\sigma_{\overline{v_{\phi,j}}}^{2}}\right)\right],$
(8)

where $j$ is iterated over $R$ bins. We note that this equation assumes that
rotational and azimuthal velocities in each of the bins are independent of
each other. The terms $\overline{v_{\phi}}_{,j}$ and
$\sigma_{\overline{v_{\phi}}_{,j}}^{2}$ are obtained from the data and are the
mean and variance of the azimuthal velocity in the j-th bin respectively.

For the prior distribution in (7), we defined flat priors, where the circular
velocities are allowed within a range of $[-400,400]$ km/s. In addition to the
velocities, we defined naive priors for the scale length parameters
$h_{R}=3\pm 1$ kpc and $h_{\sigma}=21\pm 1$ so to encompass values in the
literature (Eilers et al. 2019; Bland-Hawthorn & Gerhard 2016). As mentioned
previously, the Galactocentric distance $R_{0}$ was also treated as a free
parameter of the analysis and was given a uniform prior within $[7.8,8.5]\,\rm
kpc$. Since the solar Galactocentric velocities $V_{\odot}$ and $W_{\odot}$
are scaled with $R_{0}$ (see Eq. (1)), the chosen prior is reflected both in
the median value and error bar of the circular velocity in a given bin.

Having defined our likelihood and prior functions to use in the fitting, we
set up our MCMC algorithm using the python package emcee (Foreman-Mackey et
al. 2013). The parameter space of our model was explored using 48 independent
walkers. All in all, we used 13 parameters in the fitting, where the first ten
were circular velocities $v_{c,j}$ of the radial bins and the rest were the
Sun’s Galactocentric distance and the scale length terms.

By treating the Sun’s Galactocentric distance as a free parameter of the
analysis we were required to repeat the coordinate and velocity transformation
at each step in the MCMC. In addition to this, we also had to propagate the
covariance information of each star resulting in each step of the MCMC being
computationally expensive. In order to bring down the iteration time, we used
numpy (Harris et al. 2020) and cupy (Okuta et al. 2017), which make it
possible to implement the calculations on both CPU and GPU in an efficient
vectorised form.

The use of GPUs was particularly well motivated, since the parameter and
uncertainty propagation routines in our code consist largely of matrix
operations with relatively large arrays. GPU-accelerated computing libraries
(such as cupy) take advantage of the fact that modern GPUs have significantly
more threads than a CPU and are thus better at parallelising certain
computation routines than their CPU-counterparts. In the end, both numpy and
cupy were utilised simultaneously and the MCMC routine easily parallelised
across the available CPUs and GPU devices where the most computationally
demanding aspects of the pipeline were handled by the latter. The full data
was analysed by using two CPU cores per GPU and with a total of six GPUs the
computation time for each step was brought down to $\approx 11\,s$. This
translates into a 6-fold speed increase when compared to running the code with
just a single GPU and a 164-fold increase when running solely on CPUs. Using a
single GPU for the full dataset described in this work, quickly leads to
either out of memory issues or extremely long runtimes and thus it must be
noted that the feasibility of the analysis was heavily dependent on the
availability of multiple GPU devices and CPU cores. Our RGB sample of roughly
0.6 million RGB stars and the code used in our analysis can be found in zenodo
and https://github.com/HEP-KBFI/gaia-tools, respectively.

##  4 Results

###  4.1 Circular velocity curve

The circular velocity curve of our sample of RGB stars is summarised in table
1 and shown in figure 6. Inside each bin, we quote the median of the 1D
marginalised posterior probability distribution obtained in the MCMC fitting.
For the error bars, we quote the 16th and 84th percentile of the distribution.
We would like to highlight that, in the classical approach for calculating the
circular velocity curve, circular velocities in radial bins are calculated by
plugging values into equation (4). In the proposed new approach, we add a
simple kinematic model (given by (2)) on top of the Jeans equation, thus
transforming the problem into an inference procedure. This allows to introduce
nuisance parameters, such as the $h_{R}$, $h_{\sigma}$ and $R_{0}$, and
propagate their corresponding uncertainties (regardless of whether we have
normal or non-normal errors) into the final circular velocity curve via Bayes
theorem. The fitting procedure may, nonetheless, introduce additional
correlations between the radial bins. For this reason, we checked that, if we
fix the nuisance parameters $R_{0}$, $h_{R}$ and $h_{\sigma}$, the central
values of the circular velocities obtained with the new approach (MCMC
analysis) coincide with the values obtained by the classical or traditional
technique.

In figure 6, we compare our result to others from the literature. Our circular
velocity curve is in agreement with the one estimated in Ablimit et al. (2020)
using the 3D velocity vector method on $\sim 10^{3}$ classical Cepheids.
However, for $R>8\,\rm kpc$, we obtain larger circular velocities than those
calculated by the same authors but using the proper motions of $\sim 600$
classical Cepheids. The former and latter samples have around 370 Cepheids in
common and both results show that modelling assumptions and/or tracer samples
can induce differences in the estimated circular velocities of at least 10%.
We note that these changes are larger than the estimated uncertainties in this
work, which are in the $\lesssim 3$% level. Our error bars include statistical
uncertainties, which are negligible owing to the large data sample. They
further include uncertainties in the spatial-kinematic morphology of the
tracer stars (scale radius of the density profile $h_{R}$ and of the velocity
distribution $h_{\sigma}$) and in the Sun’s galactocentric distance. Circular
velocities show a mild sensitivity to $h_{R}$, specially for values $h_{R}\leq
2.5\,\rm kpc$ and, at least within the prior range explored in our analysis,
$h_{\sigma}$ and circular velocity central values are independent. On the
contrary, the adopted value of $R_{0}$ strongly affects the final circular
velocities and it is the main source of systematic uncertainties (from the
ones studied in this analysis).

In addition, our estimated circular velocities are also compatible with those
obtained in Eilers et al. (2019), Wang et al. (2022) and Zhou et al. (2023),
due to our large uncertainties compared to those estimated in these articles.
If we fixed the Sun’s galactocentric distance and total velocity in the
azimuthal direction to the values adopted in the former article, namely
$R_{0}=8.122\,\rm kpc$ and $V_{\odot}=245.8\,\rm km/s$, the estimated error
bars on the circular velocities are reduced and our results are incompatible
with Eilers et al. (2019) analysis at 1$\sigma$ for our fiducial distance
estimates (see section 4.4 for a comparison using the circular velocity curve
obtained with other distance estimates). This shows that $R_{0}$ is the main
source of uncertainty in the reconstruction of the circular velocity curve.
Moreover, if for the fixed $R_{0}$ case the prior range in the scale length
$h_{\sigma}$ is increased by a factor of three, the results remain unchanged.
In contrast, increasing the prior range in the scale length $h_{R}$ by the
same factor, the median values decrease by less than 2% and the size of error
bars remains roughly the same.

The circular velocity curve in Wang et al. (2022) was obtained by describing,
by means of the radial axisymmetric Jeans equation, the dynamics of all Gaia
DR3 stars within the region $160^{\rm o}<\ell<200^{\rm o}$ and $|Z|<3\,\rm
kpc$ that have measured radial velocities. All stars are thus described by the
same asymmetric drift. However, younger stars are expected to have a smaller
asymmetric drift than an older population of stars. In fact, Kawata et al.
(2019) estimated $v_{a}(R_{0})=0.28\pm 0.20\,\rm km/s$ using young classical
Cepheids, whether we obtain, as expected, the central larger value
$v_{a}(R_{0})=3\pm 7\rm\,km/s$ for older RGB stars. Figure 7 shows the
asymmetric drift as a function of Galactocentric distance for
$R_{0}=8.277\,\rm kpc$. The asymmetric drift mildly increases with distance to
the Galactic centre with a slope of $0.59\pm 0.12\,\rm km/s/kpc$.

Our estimated value of the circular velocity at the Sun’s position, namely
$233\pm 7\,\rm km/s$, is compared in table 2 with other estimates from the
literature. We found that the estimated gradient of the curve is extremely
sensitive to the radial interval included in its inference. If we remove the
first two radial bins where the circular velocities increase, the obtained
value is $-1.1\pm 0.3\,\rm kms^{-1}kpc^{-1}$, which points to a smooth
decrease of the circular velocities with Galactocentric distances. If we
rather include all radial bins, the estimated value of the slope is $0.4\pm
0.6\,\rm kms^{-1}kpc^{-1}$, which describes a flat circular velocity curve
within the uncertainties. Mróz et al. (2019) and Eilers et al. (2019) included
all radial intervals for the determination of the slope, and found decreasing
slopes that do not agree with the latter value. This may point out to the
presence of systematic biases in the distance estimations, as described in
section 4.4.

Table 1: Measured circular velocities $v_{c}$ , also plotted in figure 6. We quote the median of each $x$ bin, as the fitting is done using this adimensional variable, and the corresponding $R$ value is calculated for $R_{0}=8.277\,\rm kpc$ (GRAVITY Collaboration et al. 2022). x | R [kpc] |  $v_{c}$ [km/s] |  $\sigma_{v}^{-}$ [km/s] |  $\sigma_{v}^{+}$ [km/s] |  $v_{a}$ [km/s] |  $\sigma_{va}$ [km/s]
---|---|---|---|---|---|---
0.66 | 5.5 | 221.3 | 6.5 | 5.7 | 3.3 | 6.1
0.79 | 6.5 | 231.0 | 6.9 | 6.0 | 2.8 | 6.5
0.91 | 7.5 | 234.6 | 7.1 | 6.1 | 2.9 | 6.6
1.03 | 8.5 | 232.7 | 7.1 | 6.1 | 3.1 | 6.6
1.15 | 9.5 | 229.8 | 7.1 | 6.1 | 3.5 | 6.6
1.27 | 10.5 | 231.2 | 7.0 | 6.2 | 4.6 | 6.6
1.39 | 11.5 | 230.6 | 6.3 | 6.1 | 6.4 | 6.2
1.57 | 13.0 | 227.5 | 6.5 | 5.8 | 7.2 | 6.4
Figure 6: Circular velocity curve obtained from the MCMC. Grey dashed lines
have been plotted to indicate the position of each radial bin. In black (with
circles) we show the circular velocities as obtained in this paper where the
error bars correspond to the 16th and 84th percentile of the circular velocity
posterior distribution in a particular bin. We adopted $R_{0}=8.277\,\rm kpc$
to convert the adimensional coordinate $x$ into Galactocentric distance $R$.

Figure 7: Asymmetric drift profile. The black line shows the asymmetric drift
correction in each radial bin and the shaded region depicts its propagated
uncertainty. We note that the uncertainties are fully correlated between the
bins.

###  4.2 Smooth dark matter halo

In the region analysed in this article, mean radial velocities are $<5\%$ of
the azimuthal or rotational velocities. According to Chrobáková et al. (2020),
when the ratio of radial to azimuthal velocities is smaller than 10%, circular
velocities obtained through the radial axisymmetric and stationary Jeans
equation provide an unbiased estimator for the centrifugal velocities that
balance the averaged radial gravitational force. According to this result, our
circular velocity curve then provides an unbiased estimation of the
spherically averaged dynamical mass distribution within $14\,\rm kpc$.
Although this conclusion will be tested in an upcoming paper where we assess
the effects of modelling assumptions, such as the axial symmetry condition,
and the effect of Galactic substructure, using our current results we estimate
that the DM mass is

$\log_{10}\left[M_{\rm DM}(R<14\,{\rm kpc})/{\rm
M_{\odot}}\right]=11.2^{+2.0}_{-2.3}.$ (9)

Furthermore, we find that the local spherically averaged DM density is

$\rho_{\rm DM}(R_{0})=\left(0.41^{+0.10}_{-0.09}\right)\,{\rm
GeV/cm^{3}}=\left(0.011^{+0.003}_{-0.002}\right)\,{\rm M_{\odot}/pc^{3}}.$
(10)

These estimates were obtained by fitting the observed circular velocities to
the velocities predicted by Newtonian gravity for the baryonic components
(stellar bulge, disc and gas) and the DM halo. For each baryonic component we
adopted a set of three-dimensional density profiles, originally compiled in
Iocco et al. (2015). The stellar bulge mass is constrained by microlensing
towards the Galactic centre and the stellar disc is normalised by the stellar
surface density at the Sun’s position. We describe the DM distribution using a
generalised Navarro-Frenk-White density profile Zhao (1996). We compare
observed and predicted velocities using the Bayesian prescription presented in
Karukes et al. (2019, 2020). The estimates provided are Bayesian model
averages that include uncertainties in the Sun’s Galactocentric distance, the
three-dimensional density profile of bulge and disc stars, and the stellar
mass of the Galaxy.

We would like to highlight that our estimate of the local DM density is
compatible, within $1\sigma$ uncertainties, with recent estimates of this
quantity using the circular velocity curve method (Eilers et al., 2019; Mróz
et al., 2019; de Salas et al., 2019; Lin & Li, 2019; Karukes et al., 2019;
Sofue, 2020). In addition, it is compatible with local estimates using the
vertical Jeans equation (Salomon et al., 2020; Guo et al., 2020; Nitschai et
al., 2020) and with the most-recent local estimate using a novel machine
learning approach by Lim et al. (2023). Thus reinforcing the conclusion about
the spherical shape of the inner $\sim\leavevmode\nobreak\ 15$ kpc of the DM
halo, obtained by modelling stellar streams (Koposov et al., 2010; Bowden et
al., 2015) and the kinematics of halo stars Wegg et al. (2019).

###  4.3 The local standard of rest and the solar peculiar velocity

The total Galactocentric azimuthal velocity of the Sun can be used to derive
the solar peculiar velocity when we incorporate knowledge about the circular
velocity at its position. In particular, the total azimuthal velocity is often
decomposed as

$V_{\odot}=v_{c}(R_{0})+V_{\odot,LSR},$ (11)

where the last term is the Sun’s peculiar motion in the local standard of rest
(LSR). The treatment of the solar velocity as shown in equation (11) assumes
that the LSR moves in a circular orbit about the Galactic centre and
therefore, it coincides with the rotational standard of rest (RSR) in which
stars move on circular orbits in the azimuthally averaged gravitational
potential. However, in recent years it has been shown that the stellar disc
exhibits bulk motions at the kiloparsec scale (Bovy et al. 2015; Williams et
al. 2013; Khanna et al. 2023). In the presence of these large scale streaming
motions, the LSR, which is defined as the reference frame of a local
population of stars with zero velocity dispersion, might not coincide with the
RSR. Considering this, the total azimuthal velocity of the Sun can be
decomposed as (Drimmel & Poggio 2018)

$V_{\odot}=v_{c}(R_{0})+V_{LSR}+V_{\odot,LSR},$ (12)

where $V_{LSR}$ is the velocity of the LSR with respect to the RSR. This
difference of velocity between the LSR and RSR might account for the
discrepancy between locally-derived estimations of the Sun’s peculiar motion
(i.e. using the Strömberg relation) and globally-measured values using a
sample of tracers in a larger volume around the Sun. In fact, Bovy et al.
(2012) concluded that the LSR itself might not be on a circular orbit and it
is rotating $V_{LSR}\approx 12$ km/s faster than the actual RSR. This is in
agreement with the recently reported value in Khanna et al. (2023) of $\approx
10\,\rm km/s$. On the other hand, Bland-Hawthorn & Gerhard (2016) estimated
$V_{LSR}=0\pm$15\text{\,}\mathrm{k}\mathrm{m}\mathrm{/}\mathrm{s}$$.

Table 2: Circular velocity at the solar location $v_{c}(R_{0})$ as measured by different methods. In the last column, we quote the value of the Sun’s galactocentric distance that was adopted in each of the referenced articles.  Source |  $v_{c}(R_{0})$ [km/s] |  $R_{0}$ [kpc]
---|---|---
This work | $233\pm 7$ | 8.277
Zhou et al. (2023) | $234.04\pm 0.08{\rm(stat.)}\pm 1.36{\rm(sys.)}$ | $8.122\pm 0.031$
Kipper et al. (2021) | $228.4\pm 3.5$ | 8.3
Eilers et al. (2019) | $229.0\pm 0.2$ | $8.122\pm 0.031$
Kawata et al. (2018) | $236\pm 3$ | $8.2\pm 0.1$
Bobylev (2017) | $231\pm 6$ | 8
Huang et al. (2016) | $240\pm 6$ | 8.34
Bovy et al. (2012) | $218\pm 6$ | $8.1_{-0.1}^{+1.2}$

Assuming $V_{\odot,LSR}=$12.24$\,\rm km/s$ as measured in Schönrich et al.
(2010) from the Hipparcos Catalogue, we obtain $V_{LSR}=7\pm 7$ km/s. Our
estimate is still statistically compatible with zero streaming motion, but
nevertheless strengthens the hypothesis that a region around the Sun, with a
characteristic length scale of 1 kpc, exhibits a bulk motion in the azimuthal
direction of the order of 10 km/s.

###  4.4 Cautionary tale about distances

Our study is based on GSP-Phot distances (Andrae, R. et al., 2023), which have
been shown to systematically underestimate the distance beyond 2 kpc from the
Sun (Fouesneau, M. et al., 2023). This bias is also shown in the left panel of
Fig. 1, as the two-dimensional distribution of the estimated distance versus
inverse parallax is not symmetric with respect to the 1:1 line, as expected
from a Gaussian noise model for parallax measurements, but is more populated
to the right of this line. Fouesneau, M. et al. (2023) find that imposing a
cut-off in the quality of the parallax measurement of
$\varpi/\sigma_{\varpi}>10$ yields reliable heliocentric distances up to 10
kpc. And Andrae, R. et al. (2023) point out that a strict parallax quality
cut-off of $\varpi/\sigma_{\varpi}>20$ provides reliable distances. For our
sample of RGB stars, the first cut-off eliminates the systematic
underestimation of the GSP-Phot distances, although a slight overestimation of
distances appears (see left panel of Fig. 1). For this reason, as our fiducial
run, we adopted the last strict cut, which alleviates the mild overestimation.
However, in this section we describe how our results would change if we had
adopted the less stringent cut-off in parallax quality or rather used photogeo
distances from Bailer-Jones et al. (2021) (hereinafter referred to as BJ
distances).

A bias in distance estimates has a noticeable impact in the results of our
analysis. In particular, underestimated distances lead to an overestimation of
the circular velocity curve gradient, and thus to an underestimation of the DM
content, and vice versa in the case of overestimated distances. We note that
the bias in the estimated slope is more pronounced the greater the actual
slope. In order to assess the effect of biased distances, we performed our
MCMC analysis using four different distance estimates: BJ distances with a
cut-off of $\varpi/\sigma_{\varpi}>5$ and $\varpi/\sigma_{\varpi}>10$, and
GSP-Phot distances with $\varpi/\sigma_{\varpi}>10$ and
$\varpi/\sigma_{\varpi}>20$. For each of these distance estimates, the top
panel of figure 8 shows the resultant circular velocity curve while fixing the
Sun’s galactocentric distance and total velocity in the azimuthal direction to
the values adopted in Eilers et al. (2019), namely $R_{0}=8.122\,\rm kpc$ and
$V_{\odot}=245.8\,\rm km/s$, and leaving $h_{R}$ and $h_{\sigma}$ free. The
bottom panel of the same figure depicts circular velocities while letting
$R_{0}$ as an additional free parameter. From this figure, it is clear that
the inclusion of uncertainties in $R_{0}$ makes the four circular velocities
compatible within uncertainties. Furthermore, as we increase the cut-off in
the parallax quality from 5 to 10 for BJ distances, the declining of the curve
becomes less steep, thus increasing the DM mass of the Galaxy. On the other
hand, by increasing the cut-off from 10 to 20 for the GSP-Phot distances, we
are alleviating the mild overestimation of distances and the positive gradient
of the curves becomes shallower, reducing the DM mass. Table 3 summarises
these results.

Figure 8: Circular velocity curve for different distances estimates as explained in the main text. In the top panel, we show the results where fixing the Sun’s galactocentric distance to $R_{0}=8.122\,\rm kpc$ and leaving free the scale lengths of the radial spatial density and velocity dispersion of the tracer RGB sample. The orange band shows encompass the estimated circular velocity curve in Eilers et al. (2019) within 1$\sigma$. The results when additionally leaving $R_{0}$ as a free parameter are shown in the bottom panel.  Table 3: Summary of the results when having $R_{0}$, $h_{R}$ and $h_{\sigma}$ as free parameters and adopting different distance estimates. Namely, photogeo distances from Bailer-Jones et al. (2021) with a cut in parallax of $\varpi/\sigma_{\varpi}>5$ and $\varpi/\sigma_{\varpi}>10$, and GSP-Phot distances with quality parallax cuts of $\varpi/\sigma_{\varpi}>10$ and $\varpi/\sigma_{\varpi}>20$ (fiducial case). The second column shows the estimated slope of the circular velocity curve when a straight line is fitted to all bins (and in brackets the value obtained by eliminating the first two radial bins as the circular velocities increase within $\sim 5-7$ kpc), the local DM density and DM mass within 14 kpc are shown in the third and forth columns, respectively. We quote the velocity of the LSR in the last column. Distance Estimate | Slope (first 2 bins removed) $[\rm km/s/kpc]$ | $\rho_{0}\,\rm[GeV/cm^{3}]$ | $M_{DM}(R<14\,{\rm kpc})\,[M_{\odot}]$ | $V_{\rm LSR}\,\rm[km/s]$
---|---|---|---|---
BJ21 + $\varpi/\sigma_{\varpi}>5$ |  $-0.9\pm 0.3$ ($-1.58\pm 0.08$) | $0.37^{+0.08}_{-0.07}$ | $9.9^{+1.6}_{-1.9}$ | $7\pm 7$
BJ21 + $\varpi/\sigma_{\varpi}>10$ |  $-0.3\pm 0.3$ ($-1.0\pm 0.3$) | $0.39^{+0.09}_{-0.08}$ | $10.8^{+2.0}_{-1.7}$ | $6\pm 7$
GSP-Phot + $\varpi/\sigma_{\varpi}>10$ |  $0.6\pm 0.3$ ($0.2\pm 0.3$) | $0.43^{+0.07}_{-0.06}$ | $11.6^{+1.8}_{-2.0}$ | $6\pm 7$
GSP-Phot + $\varpi/\sigma_{\varpi}>20$ |  $0.4\pm 0.6$ ($-1.1\pm 0.3$) | $0.41^{+0.10}_{-0.09}$ | $11.2^{+2.0}_{-2.3}$ | $7\pm 7$

##  5 Summary and conclusions

We estimated the circular velocity curve from 5 kpc to 14 kpc from the
Galactic centre using $665\,660$ RGB stars that are approximately located in
one quarter of the stellar disc with 6D phase-space information as measured by
Gaia DR3, and GSP-Phot distance estimates. We determined the circular velocity
curve by describing observed rotational velocities, in adimensional radial
bins, as the difference between the circular velocity and the asymmetric
drift. The latter given by the stationary and axisymmetric radial Jeans
equations, under the further assumption of reflection symmetry above and below
the Galactic plane. In the traditional approach, one simply plugs values into
the Jeans equation. In our approach, by describing the observed rotational
velocity as the circular velocity minus the asymmetric drift, we transformed
the problem into an inference procedure. In particular, observed and model
rotational velocities were fitted using a Bayesian inference approach that
incorporates systematic and statistical uncertainties as nuisance parameters.
This allowed us to propagate into the final results uncertainties of different
nature. In particular, our relative uncertainties are $\sim 3\%$ and, apart
from the statistics, account for uncertainties in the Sun’s galactocentric
distance (which is the main source of uncertainty) and uncertainties in the
spatial-kinematic morphology of the stellar disc.

We studied the effect of biased distances on our results and showed, as
expected, that underestimated distances lead to steeper (negative) slopes and
thus to an underestimation of the dark matter content in the Galaxy. This may
explain some recent findings of significantly declining circular velocity
curves and, consequently, lower spherically averaged local DM densities than
those purely local values obtained using stars in the Solar neighbourhood.
Owing to the spherical shape of the DM halo in the inner $\sim$15 kpc of the
Galaxy, these two sets of estimates should converge.

###### Acknowledgements.

We thank the referee for her/his constructive comments and for pointing out
the biases in the distance estimates. This has undoubtedly opened up many
avenues for further studies. The authors would like to thank A. Cuoco, G.
Battaglia and E. Fernández Alvar for fuitful discussions and comments. This
work was supported by the Estonian Research Council grants PRG1006, PSG700,
PRG803, PSG864, MOBTP187, PRG780, MOBTT86 and by the European Regional
Development Fund through the CoE program grant TK133. This research 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. The authors gratefully
acknowledge the support of Nvidia, whose technology played a critical role in
the success of this research. GFT acknowledge support from the Agencia Estatal
de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) under
grant number CEX2019-000920-S and the AEI-MICIN under grant number
PID2020-118778GB-I00/10.13039/501100011033

## References

  * Ablimit et al. (2020) Ablimit, I., Zhao, G., Flynn, C., & Bird, S. A. 2020, ApJ, 895, L12
  * Abuter et al. (2018) Abuter, R. et al. 2018, Astron. Astrophys., 615, L15
  * Andrae, R. et al. (2023) Andrae, R., Fouesneau, M., Sordo, R., et al. 2023, A&A, 674, A27
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Benito et al. (2017) Benito, M., Bernal, N., Bozorgnia, N., Calore, F., & Iocco, F. 2017, J. Cosmology Astropart. Phys., 2017, 007
  * Benito et al. (2019) Benito, M., Cuoco, A., & Iocco, F. 2019, J. Cosmology Astropart. Phys., 2019, 033
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529
  * Bobylev (2017) Bobylev, V. V. 2017, Astronomy Letters, 43, 152
  * Bovy (2023) Bovy, J. 2023, Dynamics and Astrophysics of Galaxies (Princeton University Press, Princeton, NJ (in preparation))
  * Bovy et al. (2012) Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131
  * Bovy et al. (2015) Bovy, J., Bird, J. C., Pérez, A. E. G., et al. 2015, The Astrophysical Journal, 800, 83
  * Bowden et al. (2015) Bowden, A., Belokurov, V., & Evans, N. W. 2015, MNRAS, 449, 1391
  * Chrobáková et al. (2020) Chrobáková, Ž., López-Corredoira, M., Sylos Labini, F., Wang, H. F., & Nagy, R. 2020, A&A, 642, A95
  * de Salas et al. (2019) de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, J. Cosmology Astropart. Phys., 2019, 037
  * Do et al. (2019) Do, T., Hees, A., Ghez, A., et al. 2019, Science, 365, 664
  * Drimmel & Poggio (2018) Drimmel, R. & Poggio, E. 2018, Research Notes of the American Astronomical Society, 2, 210
  * Eilers et al. (2019) Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Fouesneau, M. et al. (2023) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2023, A&A, 674, A28
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Drimmel, R., Romero-Gómez, M., et al. 2023, A&A, 674, A37
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018, A&A, 616, A11
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * GRAVITY Collaboration et al. (2022) GRAVITY Collaboration, Abuter, R., Aimar, N., et al. 2022, A&A, 657, L12
  * GRAVITY Collaboration et al. (2020) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2020, A&A, 636, L5
  * GRAVITY Collaboration et al. (2019) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2019, A&A, 625, L10
  * GRAVITY Collaboration et al. (2021) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2021, A&A, 647, A59
  * Guo et al. (2020) Guo, R., Liu, C., Mao, S., et al. 2020, MNRAS, 495, 4828
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85
  * Hobbs et al. (2018) Hobbs, D., Lindegren, L., Bastian, U., et al. 2018, Gaia DR2 documentation Chapter 3: Astrometry, https://gea.esac.esa.int/archive/documentation/GDR2/
  * Huang et al. (2016) Huang, Y., Liu, X.-W., Yuan, H.-B., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 2623
  * Iocco et al. (2015) Iocco, F., Pato, M., & Bertone, G. 2015, Nature Physics, 11, 245
  * Johnson & Soderblom (1987) Johnson, D. R. H. & Soderblom, D. R. 1987, AJ, 93, 864
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864
  * Karukes et al. (2019) Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2019, J. Cosmology Astropart. Phys., 2019, 046
  * Karukes et al. (2020) Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2020, J. Cosmology Astropart. Phys., 2020, 033
  * Kawata et al. (2018) Kawata, D., Bovy, J., Matsunaga, N., & Baba, J. 2018, Monthly Notices of the Royal Astronomical Society, 482, 40
  * Kawata et al. (2019) Kawata, D., Bovy, J., Matsunaga, N., & Baba, J. 2019, MNRAS, 482, 40
  * Khanna et al. (2023) Khanna, S., Sharma, S., Bland-Hawthorn, J., & Hayden, M. 2023, Monthly Notices of the Royal Astronomical Society, 520, 5002
  * Kipper et al. (2021) Kipper, R., Tenjes, P., Tempel, E., & de Propris, R. 2021, MNRAS, 506, 5559
  * Kipper et al. (2016) Kipper, R., Tenjes, P., Tihhonova, O., Tamm, A., & Tempel, E. 2016, MNRAS, 460, 2720
  * Koposov et al. (2010) Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260
  * Kuzmin (1970) Kuzmin, G. G. 1970, arXiv e-prints, arXiv:2201.04136
  * Leung et al. (2022) Leung, H. W., Bovy, J., Mackereth, J. T., et al. 2022, Monthly Notices of the Royal Astronomical Society, 519, 948
  * Lim et al. (2023) Lim, S. H., Putney, E., Buckley, M. R., & Shih, D. 2023, arXiv e-prints, arXiv:2305.13358
  * Lin & Li (2019) Lin, H.-N. & Li, X. 2019, MNRAS, 487, 5679
  * Malkin (2013) Malkin, Z. M. 2013, Astronomy Reports, 57, 128
  * Mróz et al. (2019) Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10
  * Nitschai et al. (2020) Nitschai, M. S., Cappellari, M., & Neumayer, N. 2020, MNRAS, 494, 6001
  * Okuta et al. (2017) Okuta, R., Unno, Y., Nishino, D., Hido, S., & Loomis, C. 2017, in Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS)
  * Reid & Brunthaler (2020) Reid, M. J. & Brunthaler, A. 2020, ApJ, 892, 39
  * Rybizki et al. (2022) Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597
  * Salomon et al. (2020) Salomon, J.-B., Bienaymé, O., Reylé, C., Robin, A. C., & Famaey, B. 2020, A&A, 643, A75
  * Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829
  * Schönrich et al. (2010) Schönrich, R., Binney, J., & Dehnen, W. 2010, Monthly Notices of the Royal Astronomical Society, 403, 1829
  * Sofue (2020) Sofue, Y. 2020, Galaxies, 8, 37
  * Tempel & Tenjes (2006) Tempel, E. & Tenjes, P. 2006, MNRAS, 371, 1269
  * Thomas et al. (2019) Thomas, G. F., Laporte, C. F. P., McConnachie, A. W., et al. 2019, MNRAS, 483, 3119
  * Wang et al. (2022) Wang, H.-F., Žofia Chrobáková, López-Corredoira, M., & Labini, F. S. 2022, The Astrophysical Journal, 942, 12
  * Wegg et al. (2019) Wegg, C., Gerhard, O., & Bieth, M. 2019, MNRAS, 485, 3296
  * Williams et al. (2013) Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101
  * Zhao (1996) Zhao, H. 1996, MNRAS, 278, 488
  * Zhou et al. (2023) Zhou, Y., Li, X., Huang, Y., & Zhang, H. 2023, The Astrophysical Journal, 946, 73

##  Appendix A MCMC results

Figure 9 shows the marginalised two-dimensional and one-dimensional posterior
distributions for three different MCMC runs: first, $(R_{0},h_{R},h_{\sigma})$
are free nuisance parameters (black), second, ($h_{R},h_{\sigma})$ are not
fixed and $R_{0}$ is fixed to the value $R_{0}=8.122\,\rm kpc$ (red) and
finally, the output of the MCMC when all nuisance parameters are fixed to the
values $R_{0}=8.277\,\rm kpc$, $h_{R}=3\,\rm kpc$, $h_{\sigma}=21\,\rm kpc$.

For the first run (i.e. $R_{0}$, $h_{R}$ and $h_{\sigma}$ are free
parameters), the movement of circular velocities is driven by changes in
$R_{0}$. Since we are not considering strong priors on $R_{0}$, a strong
positive correlation between this parameter and the circular velocities is
observed. The circular velocity curve is sensitive to the Sun’s Galactocentric
distance $R_{0}$ (Benito et al. 2019), nonetheless, $R_{0}$ is constrained in
the literature much better by different types of analysis than the circular
velocities (see Malkin 2013 for a review of techniques). Therefore, we do not
aim to restrict $R_{0}$, but to assess the impact of its uncertain value on
the circular velocity curve. For this reason, we adopt as a prior a uniform
distribution that encompasses the most recent determinations of $R_{0}$ within
2$\sigma$ uncertainties. This is a conservative approach that does not favour
any particular estimate. We would like to also emphasise that the actual value
of $R_{0}$ may be troublesome. For example, the LMC is causing differences in
reflex motion in distinct parts of the Galaxy causing the definition of centre
to be vague or we do not know to what extent the different definitions of
centre (e.g. local isopotential curve, SMBH position) affect estimations via
the Jeans equations. One of the results of our analysis is that, given the
scatter in the most recent determination of $R_{0}$, this parameter represents
the main source of uncertainty in the calculation of circular velocities in
radial bins.

On the contrary, circular velocities show a mild sensitivity to $h_{R}$ and
$h_{\sigma}$. Neither of these scale lengths can be constrained by our
analysis and are simply treated as nuisance parameters whose prior range is
defined by observational determinations (see e.g. Eilers et al. 2019 and
references therein).

Figure 9: One and two-dimensional marginalised posterior distributions of the
circular velocity. The figure shows distributions from three different
parameter setups: all nuissance parameters free (black), only scale parameters
free (red), all nuissance parameters fixed (blue), see text for more details.
The contours delimit regions of 2-$\sigma$ probability and the best fit
circular velocity of the black posteriors is demarked with horisontal and
vertical lines. The diagonal contains the normalised 1D posteriors of the
circular velocities from the three different runs.