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

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. 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. 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\pm7\, \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.


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 in-stead 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 the-Article number, page 1 of 12 arXiv:2309.02895v1 [astro-ph.GA] 6 Sep 2023 A&A proofs: manuscript no. main orem. 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.

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 K and 5500 K and surface gravity satisfying the condition log g < 3. Both stellar parameters are provided as data products in Gaia DR3 . 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 ≤ 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 kpc (which removes ca. 4.5M stars), a cut in Galactocentric distances 5 kpc < R < 14 kpc (removing ca. 84k stars), and a cut in heliocentric velocity |v − v ⊙ | < 210 km/s. 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).
We used GSP-Phot distances  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 ϖ/σ ϖ > 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 ϖ/σ ϖ ≲ 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.

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 ⊙ , V ⊙ , W ⊙ ) 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 ∈ [7.8 − 8.5] 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, 2020Leung 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σ 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) 2 . 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 . 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 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σ level (GRAVITY Collaboration et al. 2022). In this way, we obtain the following vector where U ⊙ , V ⊙ , W ⊙ correspond to the velocity components in the Galactocentric x, y, and z-directions, respectively. As U ⊙ 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 ϖ/σ ϖ > 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, ϕ, z, v R , v ϕ , v z ).  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 kpc, 12 kpc and 14 kpc to ease visualisation. For the transformation to Galactocentric coordinades, we adopt R 0 = 8.277 kpc, z 0 = 25 pc and v ⊙ = (11.1, 251.5, 8.59) km/s.
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.

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 ). Furthermore, we note that by marginalising over the azimuthal coordinate ϕ, 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 ∆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.

Axisymmetric kinematic model
Inside each radial bin, we modeled the mean rotational or azimuthal velocity as 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 Substituting v 2 ϕ = σ * 2 ϕ + (v ϕ ) 2 iii and assuming v 2 R = σ * 2 R , we further obtain iii To 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. We have checked the latter assumption explicitly and observed that the measured values of v 2 R are at least two orders of magnitude larger than the measured 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 The sum v c + v ϕ in the denominator is often approximated as ≈ 2v c (Binney & Tremaine 2008). Nonetheless, we leave it as it is and estimate v ϕ as the mean rotational velocity inside each bin. In addition to this, the diagonal components of the velocity-dispersion tensor σ * 2 ϕ and σ * 2 R 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 ν is described by an exponential profile, namely ν ∝ exp (−R/h R ) with h R the disc scale radius. Notice that in the 4th component we describe σ * R as σ * R ∝ exp(−R/h σ ), where h σ 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 We neglect the last term of (5) in our axisymmetric treatment of the rotation curve derivation as v R v z ≈ 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 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).

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 θ, namely the circular velocities {v c, j }, h R , h Θ and R 0 .
According to Bayes' theorem, the posterior distribution of a set of model parameters θ given a particular set of data D can be defined as where p(θ) 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 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 v ϕ , j and σ 2 v ϕ , j 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 ± 1 kpc and h σ = 21 ± 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] kpc. Since the solar Galactocentric velocities V ⊙ and W ⊙ 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 ≈ 11 s. This translates into a 6-fold speed increase when compared to running the code with just a single GPU and a 164fold 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.

Circular velocity curve
The circular velocity curve of our sample of RGB stars is summarised in table 4.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 σ 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 σ , 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 ∼ 10 3 classical Cepheids. However, for R > 8 kpc, we obtain larger circular velocities than those calculated by the same authors but using the proper motions of ∼ 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 ≲ 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 σ ) and in the Sun's galactocentric distance. Circular velocities show a mild sensitivity to h R , specially for values h R ≤ 2.5 kpc and, at least within the prior range explored in our analysis, h σ 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 kpc and V ⊙ = 245.8 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σ for our fiducial distance esti-mates (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 σ 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 o < ℓ < 200 o and |Z| < 3 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 ± 0.20 km/s using young classical Cepheids, whether we obtain, as expected, the central larger value v a (R 0 ) = 3 ± 7 km/s for older RGB stars. Figure 7 shows the asymmetric drift as a function of Galactocentric distance for R 0 = 8.277 kpc. The asymmetric drift mildly increases with distance to the Galactic centre with a slope of 0.59 ± 0.12 km/s/kpc.
Our estimated value of the circular velocity at the Sun's position, namely 233 ± 7 km/s, is compared in table 4.3 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 ± 0.3 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 ± 0.6 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.

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 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 M DM (R < 14 kpc)/M ⊙ = 11.2 +2.0 −2.3 .
Furthermore, we find that the local spherically averaged DM density is ρ DM (R 0 ) = 0.41 +0.10 −0.09 GeV/cm 3 = 0.011 +0.003 −0.002 M ⊙ /pc 3 . (10) These estimates were obtained by fitting the observed circular velocities to the velocities predicted by Newtonian gravity for 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 kpc (GRAVITY Collaboration et al. 2022   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. (2019Karukes et al. ( , 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σ 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 ∼ 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).

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 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) where V LS R 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 LS R ≈ 12 km/s faster than the actual RSR. This is in agreement with the recently reported value in Khanna et al. (2023) of ≈ 10 km/s. On the other hand, Bland-Hawthorn & Gerhard (2016) estimated V LS R = 0 ± 15 km/s. Assuming V ⊙,LS R =12.24 km/s as measured in  from the Hipparcos Catalogue, we obtain V LS R = 7 ± 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.

Cautionary tale about distances
Our study is based on GSP-Phot distances , which have been shown to systematically underestimate the distance beyond 2 kpc from the Sun ). 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 ϖ/σ ϖ > 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 ϖ/σ ϖ > 20 provides reliable distances. For our sample of RGB stars, the first cutoff 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 ϖ/σ ϖ > 5 and ϖ/σ ϖ > 10, and GSP-Phot distances with ϖ/σ ϖ > 10 and ϖ/σ ϖ > 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 kpc and V ⊙ = 245.8 km/s, and leaving h R and h σ 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 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.  Table 3. Summary of the results when having R 0 , h R and h σ as free parameters and adopting different distance estimates. Namely, photogeo distances from Bailer-Jones et al. (2021) with a cut in parallax of ϖ/σ ϖ > 5 and ϖ/σ ϖ > 10, and GSP-Phot distances with quality parallax cuts of ϖ/σ ϖ > 10 and ϖ/σ ϖ > 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 ∼ 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.

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 ∼ 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 So-lar neighbourhood. Owing to the spherical shape of the DM halo in the inner ∼15 kpc of the Galaxy, these two sets of estimates should converge.