# Detection of the Keplerian decline in the Milky Way rotation curve

Yongjun Jiao, Francois Hammer, Haifeng Wang, Jianling Wang, Philippe

Amram, Laurent Chemin, Yanbin Yang

► To cite this version:

Yongjun Jiao, Francois Hammer, Haifeng Wang, Jianling Wang, Philippe Amram, et al. Detection of the Keplerian decline in the Milky Way rotation curve. Astronomy & Astrophysics - A&A, 2023, 678, pp.A208. 10.1051/0004-6361/202347513. hal-04211216

HAL Id: hal-04211216

https://hal.science/hal-04211216v1

Submitted on 9 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Detection of the Keplerian decline in the Milky Way rotation curve

Yongjun Jiao<fnref n="1" />, François Hammer<fnref n="1" />, Haifeng Wang<fnref n="2" />, Jianling Wang<fnref n="1" />,<fnref n="3" />, Philippe Amram<fnref n="4" />,
Laurent Chemin<fnref n="5" />, and Yanbin Yang<fnref n="1" />

Received 20 July 2023 / Accepted 21 August 2023

## ABSTRACT

Our position inside the Galactic disc has previously prevented us from establishing an accurate rotation curve (RC). The advent of Gaia and its third data release (Gaia DR3) made it possible to specify the RC up to twice the optical radius. We aim to establish a new RC of the Galaxy from the Gaia DR3 by drastically reducing systematic uncertainties. Our goal is to provide a new estimate of the mass of the Galaxy. We compared different estimates, established a robust assessment of the systematic uncertainties, and addressed differences in methodologies, particularly regarding distance estimates. We find a sharply decreasing RC for the Milky Way; the decrease in velocity between 19.5 and 26.5 kpc is approximately 30 km s⁻<fnref n="1" />. We identify, for the first time, a Keplerian decline of the RC, starting at ~19 kpc and ending at ~26.5 kpc from the Galaxy centre, while a flat RC is rejected with a significance of 3σ. The total mass is revised downwards to 2.06±0.13 × 10¹<fnref n="1" /> M☉, which is in agreement with the absence of a significant mass increase at radii larger than 19 kpc. We evaluated the upper limit on the total mass by considering the upper values of velocity measurements, which leads to a strict, unsurpassable limit of 5.4 × 10¹<fnref n="1" /> M☉.

**Key words.** Galaxy: kinematics and dynamics – Galaxy: general – Galaxy: stellar content – Galaxy: structure

## 1. Introduction

Almost a century ago, Lundmark (1925, see their Fig. 4) was amongst the first to identify the flat behaviour of disc-galaxy rotation curves (RCs). From optical spectroscopy, Babcock (1939, then Mayall 1951) reported that the RC of M31 shows no decrease up to 20 kpc from the centre. With a larger sample of galaxies, Rubin et al. (1978) found that several spiral galaxies have a flat RC. The advent of radio astronomy made it possible to probe galaxy rotation beyond the optical disc. The increased sensitivity of radio telescopes allowed Bosma (1978) to obtain the first sample of galaxies observed in the neutral hydrogen line, and he demonstrated that most galaxies show a flat RC. Extended flat RCs of spiral galaxies can be considered as major evidence of the presence of an extended halo of dark matter (DM) surrounding them.

Early RC investigations were made beyond our own Galaxy because our location in it and the extinction prevented straightforward, direct determination of the RC as is possible for nearby external galaxies. HII regions, OB stars, carbon stars, planetary nebulae, and cepheids have been used to trace the rotation of the Galaxy. As reviewed by Schmidt (1965), before the early 1960s, the outer RC was thought to be Keplerian, and from the late 1960s until today, evidence has been found to show that its outer curve is rather flat, albeit with significant uncertainties. Therefore, the Galaxy must contain large amounts of dark matter. However, it is interesting to note that the RC collected by Sofue et al. (2009) is, despite huge uncertainties, consistent with a decreasing RC from ~15 to ~23 kpc.

A revolution in this area came with Gaia, whose proper motion measurements have allowed 3D velocity measurements. Pioneering efforts to establish the Milky Way (MW) RC with Gaia DR2 were made by Eilers et al. (2019), who used spectrophotometric parallax distances from Hogg et al. (2019) for a set of more than 20 000 red giant branch (RGB) stars. Eilers et al. (2019) provided by far the most accurate RC for the MW, a success confirmed by Mróz et al. (2019) on the basis of a small sample of variable stars that nevertheless led to very accurate distance estimates.

The MW has been found to have an exceptionally quiet merger history, evidenced by its rather pristine halo and small angular momentum when compared to other spiral galaxies (Hammer et al. 2007, and references therein). This was confirmed by the discovery of its last major merger, Gaia-Sausage-Enceladus (GSE, 8–10 Gyr ago, Belokurov et al. 2018; Haywood et al. 2018; Helmi et al. 2018), which was identified from the angular momentum signatures of residual stars in the halo. Because most spiral galaxies underwent their last major merger more recently (~6 Gyr ago, Hammer et al. 2009; Puech et al. 2012), the MW disc may be less affected by the large-scale motions expected following such a major event. This would make the MW one of the most appropriate targets for deriving a RC, assuming equilibrium conditions. It has been argued that a minor merger, such as the infall of Sagittarius (Sgr) 4–6 Gyr ago, could have perturbed the MW outer disc, providing a possible explanation for the warp (Bailin 2003) and the observed the vertical oscillations (Laporte et al. 2018). Being at first passage, the Large Magellanic Cloud

<footnote n="1">¹ GEPI, Observatoire de Paris, Paris Sciences et Lettres, CNRS, Place Jules Janssen, 92195 Meudon, France
e-mail: yongjun.jiao@obspm.fr, francois.hammer@obspm.fr</footnote>

<footnote n="2">² Centro Ricerche Enrico Fermi, Via Panisperna 89a, 00184 Rome, Italy</footnote>

<footnote n="3">³ CAS Key Laboratory of Optical Astronomy, National Astronomical Observatories, Beijing 100101, PR China</footnote>

<footnote n="4">⁴ Aix-Marseille Univ., CNRS, CNES, LAM, 38 Rue Frédéric Joliot Curie, 13338 Marseille, France</footnote>

<footnote n="5">⁵ Instituto de Astrofísica, Universidad Andes Bello, Fernandez Concha 700, Las Condes, Santiago, RM, Chile</footnote>
(LMC) may have affected the MW disc location with respect to its halo (Conroy et al. 2021; Erkal et al. 2021), but is unlikely to have affected its internal dynamics or morphology (Laporte et al. 2018).

Gaia DR3 (Gaia Collaboration 2023b) has provided improved parallaxes and proper motions, whose systematic uncertainties are smaller by a factor of 2 when compared to those of Gaia DR2 (Lindegren et al. 2021a,b). Gaia Collaboration (2023a, hereafter D23) measured a robust RC of the MW out to R = 14 kpc from the 3D velocity space of a clean sample of OB stars and another sample of RGB stars using Gaia DR3 data. This was followed by three studies investigating the outer MW RC (Wang et al. 2023a; Zhou et al. 2023; Ou et al. 2023), which used different star samples and methodologies. Wang et al. (2023a) derived distances from the Gaia parallaxes (π) of a very large number of stars for which there are also radial velocities from Gaia DR3. These latter authors implemented a Lucy's inversion method (LIM; Lucy 1974), from which they populated 6D phase space cells (l, b, π, V<sub>r</sub>, μ<sub>α</sub><sup>*</sup>, μ<sub>δ</sub>). These authors then determined the average values of the velocity components and their dispersions, which were tested and shown to provide reliable results. Ou et al. (2023) followed a technique similar to that of Eilers et al. (2019), namely using spectrophotometric distances to establish the RC of over 30 000 RGB stars. A similar approach was adopted by Zhou et al. (2023), although their distances were estimated using different priors.

Figure 1 compares the three corresponding RCs; they show many consistencies, except in the inner (where Wang et al. 2023a have lower velocities) and outer (where Zhou et al. 2023 have higher velocities) parts. This latter discrepancy should significantly affect the estimation of the dynamical mass<fnref n="1" /> of the MW. The goal of the present paper is to verify whether these differences can be attributed to different methodologies in order to establish the MW RC and to determine the dynamical mass of our Galaxy from Gaia DR3. In Sect. 2, we describe the Jeans equation, and show how we derive the systematic uncertainties for the Wang et al. (2023a) RC using a method similar to that of Eilers et al. (2019) and to that of Ou et al. (2023). In Sect. 3, we show the approach we take to determine the MW dynamical mass using a set of models including baryonic and DM components. In Sect. 4, we show that the three RCs can be reconciled together, providing a sharp velocity drop at large radii, and we provide the MW dynamical mass range and its uncertainty. In Sect. 5, we describe the limitations of our estimates, and compare them to estimates based on other tracers. We also test whether or not we can detect the Keplerian decline of the MW RC.

## 2. Methods: RC determination and associated uncertainties

### 2.1. The RC derived from Jeans equation

Assuming an axisymmetric MW potential and a disc at equilibrium, Wang et al. (2023a) used the Jeans equation (Binney & Tremaine 2008, Eq. (4.222a)) to measure the circular velocity curve:

$$ \frac{\partial v \langle V_R \rangle}{\partial t} + \frac{\partial v \langle V_R^2 \rangle}{\partial R} + \frac{\partial v \langle V_R V_z \rangle}{\partial z} + v \left( \frac{\langle V_R^2 \rangle - \langle V_\phi^2 \rangle}{R} + \frac{\partial \Phi}{\partial R} \right) = 0, \quad (1) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of the three different measurements of the MW RC based on Gaia DR3.</figcaption></figure>

where v denotes the matter density distribution. By assuming a steady state, a disc that is symmetric about its equator (∂v/∂z = 0 at z = 0), and an exponential radial profile of the tracer population with a scale length of h<sub>R</sub> = 2.5 kpc (Jurić et al. 2008), the circular velocity curve can be derived using R(∂Φ/∂R) = V<sub>c</sub><sup>2</sup> in Eq. (1):

$$ V_c^2 = \langle V_\phi^2 \rangle + \frac{R - h_R}{h_R} \langle V_R^2 \rangle - R \frac{\partial \langle V_R^2 \rangle}{\partial R} - R \frac{\partial \langle V_R V_z \rangle}{\partial z}. \quad (2) $$

For the three velocity components, we apply:

$$ \langle V_X^2 \rangle = \langle V_X \rangle^2 + \sigma_{\langle V_X \rangle}^2, \quad (3) $$

where X presents R, φ, and z, respectively.

Concerning the RC of Wang et al. (2023a), we adjusted the bin size to improve the calculation within each bin. Specifically, we calculated ⟨V<sub>φ</sub><sup>2</sup>⟩ and ⟨V<sub>R</sub><sup>2</sup>⟩ using the same 1 kpc width bins as was done by Wang et al. (2023a). However, we derived the gradient term ∂⟨V<sub>R</sub><sup>2</sup>⟩/∂R after centring the 1 kpc bin width around each data point. Figure 2 shows that the slightly modified RC is consistent with that of Wang et al. (2023a) within 22 kpc, which means that systematic uncertainties associated with the choice of the bin size and position are very small. The largest modification is for the point at R = 24.5 kpc, whose velocity amplitude decreases by 8 km s<sup>-1</sup>. However, this velocity is more consistent with the decreasing slope of the RC between 13 and 23 kpc. Data beyond 27.5 kpc are lacking (within a height of 3 kpc, i.e. |z| < 3 kpc), and so we do not derive the last point at R = 27.5 kpc.

We derived the RC and the associated systematic uncertainties within the height of 3 kpc, i.e., |z| < 3 kpc, which provides the most extended RC in Wang et al. (2023a). The cross term ∂⟨V<sub>R</sub>V<sub>z</sub>⟩/∂z was also neglected in the calculations of several previous works (e.g. Eilers et al. 2019; Wang et al. 2023a) because it is generally considered to be two to three orders of magnitude smaller than the remaining terms at many radii. This term is further analysed in Sect. 2.2.1.

### 2.2. Systematic uncertainties

The MW has various and complex dynamical structures, which implies that the assumptions of the time-independent gravitational potential and of a smooth density distribution are rough approximations of the true dynamics. For example, the observed velocity field in our Galaxy has been found to have asymmetrical motions with significant gradients in all velocity components

<footnote n="1">¹ In the following, we use the term dynamical (or total) mass, (M<sub>dyn</sub>), to refer to the sum of the baryonic and dark matter masses. This latter is also referred to as the virial dark-matter mass (M<sub>vir</sub>) inside the virial radius, above which no additional mass is expected.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison of the RC of this work with that of Wang et al. (2023a). Red points and error bars represent the RC measurement with adjusted bin sizes and systematic uncertainties at a vertical height of $|z| < 3$ kpc. Grey points and error bars present previous findings from Wang et al. (2023a) without systematic uncertainties.</figcaption></figure>

(D23) because of various gravitational influences, such as those of the bar, bulge, and spiral arms, or because of the tidal interaction with the Sagittarius dwarf galaxy (Bailin 2003). It is crucial to estimate the systematic uncertainties brought to the circular velocity curve from the data. Wang et al. (2023a) derived the kinematic maps of different velocity components from Gaia DR3 up to 30 kpc. However, the perturbations due to the radial velocity component are sufficiently small, which justifies the use of the time-independent Jeans equation.

### 2.2.1. The neglected cross-term

The vertical gradient of the cross-term $\langle V_R V_z \rangle$ is usually considered negligible within about 20 kpc. As we consider the RC out to 27 kpc, the largest contribution to the systematic uncertainty of the circular velocity comes from the neglected cross-term in the Jeans equation at large radii. Figure 3 shows the map of this cross-term on the projection of the Galactic (R, z) plane, which suggests significant variations of the cross-term with radius. The red dotted curve of Fig. 4 presents the corresponding contribution to systematic uncertainties. This term causes systematic uncertainties of smaller than 2% below 23 kpc, and of up to ~8% at 25.5 kpc. It should not be neglected at large radii. In addition, we also limited the analysis to a narrower Galactic plane, that is $|z| < 2$ kpc. In this case, systematic uncertainties from this term are smaller than 5% within 25 kpc.

### 2.2.2. Disc scale length

Another contribution to the systematic uncertainty that cannot be ignored is the unknown density profile of the tracer population, particularly in the outer disc. Our calculation assumes an exponential density profile with a scale length of $h_R = 2.5$ kpc. Following Eilers et al. (2019), we vary this scale length, that is, $\Delta h_R = \pm 1$ kpc, and this causes systematic uncertainties at the ~1% level, which are represented in Fig. 4 by an orange-dashed line.

### 2.2.3. Disc radial-density profile

The functional form chosen for the density profile can also lead to systematic uncertainties. We therefore applied a power-law density profile instead, for which we chose an index of $\alpha = -2.25$, which has the same slope as the exponential function

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Cross-term $\langle V_R V_z \rangle$ using LIM in the Galactic (R, z) plane.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Summary of potential systematic uncertainties in the circular velocity curve at a vertical height of $|z| < 3$ kpc. We estimate systematic uncertainties arising from the neglected cross-term (see the red-dotted line), from varying the exponential scale length of the density profile (orange-dashed line), from passing from an exponential to a power law (blue-dashed line) with an index $\alpha = -2.25$, and finally after splitting the sample into two parts (green-dotted line). The total systematic uncertainty (black-solid line) is at the 2% level up to $R = 22.5$ kpc, and then increases at larger radii.</figcaption></figure>

at the Sun's location. This leads to a systematic error of less than 4% (blue-dashed line in Fig. 4).

### 2.2.4. Splitting data sample

In order to estimate global systematic uncertainties arising from the data sample of different azimuthal ranges, we divided the data sample into two subsamples selected within $160^{\circ} \le l \le 180^{\circ}$ and $180^{\circ} \le l \le 200^{\circ}$, respectively. We calculated the average of the velocity difference between these two subsamples and the total sample. We present this systematic contribution in Fig. 4 with a green-dotted line. The systematic uncertainty due to the data sample is smaller than 2% within 22 kpc. At large radii, it is comparable to the systematic uncertainties associated with the neglected cross-term.

### 2.2.5. Total uncertainties including systematics

The total relative systematic uncertainties are indicated by the black-solid line of Fig. 4, which is obtained from the quadratic
<table><caption>**Table 1.** Parameters of bulge, see more details in Sect. 3.1.</caption><thead><tr><th>Parameter</th><th>Bulge E</th><th>Bulge G</th><th>Bulge B2</th></tr></thead><tbody><tr><td>M<sub>bulge</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>1.962</td><td>1.639</td><td>1.550</td></tr><tr><td>x<sub>b</sub> (kpc)</td><td>0.899</td><td>1.239</td><td>0.700<sup>(a)</sup></td></tr><tr><td>y<sub>b</sub> (kpc)</td><td>0.386</td><td>0.609</td><td>—</td></tr><tr><td>z<sub>b</sub> (kpc)</td><td>0.250</td><td>0.438</td><td>—</td></tr></tbody></table>

Notes. (a) Value of $r_b$ in Eq. (7).

<table><caption>**Table 2.** Parameters of disk, see more details in Sect. 3.1.</caption><thead><tr><th>Parameter</th><th>Disk CM</th><th>Disk J</th><th>Disk dJ</th><th>Disk B2</th></tr></thead><tbody><tr><td>M<sub>thin</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>3.11</td><td>3.17</td><td>3.33</td><td>3.65</td></tr><tr><td>M<sub>thick</sub> (10<sup>10</sup> M<sub>☉</sub>)</td><td>0.82</td><td>0.90</td><td>0.78</td><td>—</td></tr><tr><td>L<sub>thin</sub> (kpc)</td><td>2.75</td><td>2.60</td><td>2.60</td><td>2.35</td></tr><tr><td>L<sub>thick</sub> (kpc)</td><td>4.10</td><td>3.60</td><td>4.10</td><td>—</td></tr><tr><td>H<sub>thin</sub> (kpc)</td><td>0.25</td><td>0.30</td><td>0.25</td><td>0.14</td></tr><tr><td>H<sub>thick</sub> (kpc)</td><td>0.75</td><td>0.90</td><td>0.75</td><td>—</td></tr></tbody></table>

sum of all systematic uncertainties assuming that the errors are all Gaussian:

$$ \sigma_{\text{sys}}^2 = \sigma_{\text{Cross Term}}^2 + \sigma_{\text{Scale Length}}^2 + \sigma_{\text{Density Profile}}^2 + \sigma_{\text{Splitting Sample}}^2 \quad (4) $$

This total has been adopted to calculate error bars in Fig. 2. The latter (see red vertical error bars) corresponds to the quadratic sum of statistical errors with systematic uncertainties.

However, we note that the systematic uncertainty due to the neglected cross-term and that obtained after splitting the sample into two parts follow the same trend (compare the red-dotted with the green-dashed line in Fig. 4). A quite similar trend may be identified in Fig. 5 of Ou et al. (2023). This leads us to suspect that systematic uncertainties are overestimated by both studies.

# 3. Methods: MW mass models

## 3.1. Varying baryonic models for bulge and disc

The contribution of the baryonic components to the MW mass is still uncertain, and this may affect the determination of the DM distribution (Karukes et al. 2020; Jiao et al. 2021). The basic idea is to cope with uncertainties in baryons by using a very large grid of possible models. However, we are aware that some baryonic models may not be fully consistent with other important constraints from the vertical dynamics of the disc stars (Bovy & Rix 2013) or from microlensing (Wegg et al. 2016).

To verify how the baryonic mass distribution may affect the RC, we followed Jiao et al. (2021) in considering several possible combinations of models for the bulge and the disc presented by Iocco et al. (2015), all with baryonic mass smaller than $7 \times 10^{11} M_{\odot}^2$. In addition, we also considered the baryonic model B2 from de Salas et al. (2019), which includes neutral gas and dust components, and was used by Ou et al. (2023). Two triaxial density profiles E (Exponential-type) and G (Gaussian-type) for the bulge from Stanek et al. (1997) can be expressed as:

$$ E: \rho_{\text{bulge}}(x, y, z) = \rho_0 e^{-r} \qquad (5) $$

$$ G: \rho_{\text{bulge}}(x, y, z) = \rho_0 e^{-r^2/2} $$

with:

$$ r^2 = \frac{x^2}{x_b^2} + \frac{y^2}{y_b^2} + \frac{z^2}{z_b^2}, \qquad (6) $$

where $(x, y, z)$ are the coordinates along the major, intermediate, and minor axes. The bulge of B2 is modelled with a Hernquist potential:

$$ \Phi(r) = -\frac{GM}{r + r_b}. \qquad (7) $$

For the disc component, we adopted a double exponential as described below:

$$ \rho(R, z) = \rho_0 \exp\left(-\frac{R}{L} - \frac{z}{H}\right), \qquad (8) $$

where $\rho_0 = M/(4\pi HL^2)$ is the normalisation, M is the corresponding disc mass, and L and H are the disc scale length and height, respectively. Three-disc models (CM from Calchi Novati & Mancini 2011, dJ from de Jong et al. 2010 and J from Jurić et al. 2008) contain thin and thick discs. Model B2 possesses only a thin disc. The dust component and the gas distribution of B2 are also modelled as double exponential profiles (Eq. (8)).

Tables 1 and 2 provide the various parameters adopted to represent the bulge and the disc of the MW.

## 3.2. Dark matter model for the halo

The Einasto profile (Einasto 1965; Retana-Montenegro et al. 2012) is widely used to describe the distributions of stellar light and of the mass of dark matter in galaxies, whose density is defined as:

$$ \rho(r) = \rho_0 \exp\left[-\left(\frac{r}{h}\right)^{1/n}\right], \qquad (9) $$

where n is the Einasto index, which can determine how fast the density decreases with r. Several studies (e.g. Chemin et al. 2011; Jiao et al. 2021) have shown that the Einasto profile gives a significantly better fit to the RCs when compared to the Navarro, Frenk & White profile (NFW; Navarro et al. 1997). In particular, Jiao et al. (2021) found that the three-parameter Einasto profile may account for a larger range of outer slopes and generate a plausible wide mass range. This is also supported by Ou et al. (2023) and Sylos Labini et al. (2023), who found the NFW profile unsuitable for reproducing the external slope of the MW RC. Therefore, the present study focuses on a spherical Einasto profile for modelling the MW RC.

Constraints on the shape of the dark matter halo in the Milky Way are weak and no consistent picture has yet emerged. However, at 20 kpc scales, Küpper et al. (2015) and Koposov et al. (2010) determined the flattening of the dark halo to be $q_z = 0.95 \pm 0.15$, that is, almost spherical. In the inner regions of the MW, the disc and bulge dominate the RC in that they contribute ~85% of the rotational velocity and ~70% of the rotational support at 2.2 disc scale lengths, i.e., where the disc velocities are maximal and therefore the disc is maximal (Sackett 1997). The wide variety of baryonic models used in this study allows us to test the hypothesis of the maximum disc and to resolve the debate over the core or cusp nature of the dark halo of our galaxy. It also allows us to test the NFW halo, which in any case cannot fit a decreasing RC.

The optical diameter of an external galaxy is often defined using the $D_{25}$ isophote. This limit has been used to estimate the radius of the MW to 13.4 kpc (e.g. Goodwin et al. 1998). In the following, we obtain a radius that is almost twice as large.

<footnote n="2">² We notice that Model I of Pouliasis et al. (2017) used by Eilers et al. (2019) and Jiao et al. (2021) overestimates the bulge mass (Paola di Matteo, priv. comm.).</footnote>
<table><caption>**Table 3.** Measurements of the Circular Velocity of the Milky Way.</caption><thead><tr><th>R [kpc]</th><th>V<sub>c</sub> [km s<sup>-1</sup>]</th><th>σ<sub>V<sub>c</sub></sub> [km s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>9.5</td><td>221.75</td><td>3.17</td></tr><tr><td>10.5</td><td>223.32</td><td>3.02</td></tr><tr><td>11.5</td><td>220.72</td><td>3.47</td></tr><tr><td>12.5</td><td>222.92</td><td>3.19</td></tr><tr><td>13.5</td><td>224.16</td><td>3.48</td></tr><tr><td>14.5</td><td>221.60</td><td>4.20</td></tr><tr><td>15.5</td><td>218.79</td><td>4.75</td></tr><tr><td>16.5</td><td>216.38</td><td>4.96</td></tr><tr><td>17.5</td><td>213.48</td><td>6.13</td></tr><tr><td>18.5</td><td>209.17</td><td>4.42</td></tr><tr><td>19.5</td><td>206.25</td><td>4.63</td></tr><tr><td>20.5</td><td>202.54</td><td>4.40</td></tr><tr><td>21.5</td><td>197.56</td><td>4.62</td></tr><tr><td>22.5</td><td>197.00</td><td>3.81</td></tr><tr><td>23.5</td><td>191.62</td><td>12.95</td></tr><tr><td>24.5</td><td>187.12</td><td>8.06</td></tr><tr><td>25.5</td><td>181.44</td><td>19.58</td></tr><tr><td>26.5</td><td>175.68</td><td>24.68</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Circular velocity of the Milky Way. The red data points are the measurements computed in this work; error bars include systematic uncertainties. The black solid line represents the sum of the baryonic and dark matter components: the baryonic model B2 (blue-dashed line), including its decomposition into baryonic components (bulge, disc, gas, and dust) and the best fit of the Einasto dark matter profile (red-dashed line).</figcaption></figure>

# 4. Results

## 4.1. Measurement of the rotation curve and comparison with Ou et al. (2023)

We report our final circular velocity in Table 3. Figure 5 shows our final circular velocity curve and one of our best fits with $n = 0.43$, $h = 11.41$ kpc, and $\rho_0 = 0.01992 M_\odot \text{ pc}^{-3}$ (see first line of Table 4). In our study, as well as in Ou et al. (2023), error bars account for systematic uncertainties. For $R \ge 13$ kpc, the two RCs are in reasonably good agreement except for one point at about 23 kpc (see Fig. 1). We suspect that this discrepancy is caused by the disagreement over the radial velocity component at $R \sim 23$ kpc, for which the top panel of Fig. 3 of Ou et al. (2023) shows a large deviation on $\sqrt{\langle v_R^2 \rangle}$. In the range of $R = 9-13$ kpc, our RC points are slightly lower than those of Ou et al. (2023), which is discussed in Appendix A.

The largest discrepancy between the RC of this paper and that of Ou et al. (2023) is perhaps related to the amplitude of the error bars, which are larger in this latter study (compare Fig. 4 with Fig. 5 of Ou et al. 2023).

Both RCs show a significant decline with increasing radius, which can be well approximated by a linear function (see Fig. 2):

$$V_c(R) = V(R_\odot) + \beta(R - R_\odot), \quad (10)$$

where $R_\odot$ is the distance between the Sun and the Galactic centre<fnref n="3" />. We find that the slope of our declining RC is $\beta = -(2.18 \pm 0.23) \text{ km s}^{-1} \text{ kpc}^{-1}$, which is similar to the value of $\beta = -(2.22 \pm 0.20) \text{ km s}^{-1} \text{ kpc}^{-1}$ obtained by Ou et al. (2023)<fnref n="4" />.

Wang et al. (2023a) also split the Galactic region into two, one with galactic latitude $b > 0^\circ$ and the other with $b < 0^\circ$ (or

one with $z > 0$ kpc and the other with $z < 0$ kpc) and found an uncertainty on the slope of RC of ~20%.

## 4.2. Comparison with Zhou et al. (2023)

Figure 1 shows that the RC from Zhou et al. (2023) indicates larger velocities at the MW disc outskirts. In Appendix B, we compare the distances adopted by the present study to those adopted in other estimates (see Fig. B.1), which leads us to suspect that the distances by Zhou et al. (2023) are overestimated. After correcting for this, it appears that the Zhou et al. (2023) RC is consistent with both the RC of the present study and that from Ou et al. (2023). We also notice that Zhou et al. (2023) did not consider the impact of the cross-term when analysing systematic uncertainties. For consistency, we have not considered this study in the following.

## 4.3. Estimated range for the dynamical mass of the Milky Way

Using a Bayesian analysis, one can determine the posterior distribution of the model parameters based on the given data. In the present study, we applied the Markov chain Monte Carlo (MCMC) affine invariant sampler EMCEE<fnref n="5" /> (Foreman-Mackey et al. 2013) to test the parameter space of the Einasto profile using flat priors; that is, $M_0 = 4\pi h^3 \rho_0$, $h$, and $1/n$, from $10^{10}$ to $10^{14} M_\odot$, from 0 to 20, and from 0 to 5, respectively. Following previous studies, the sum of the logarithm of the likelihood for the observed RC can be derived as:

$$\ln \mathcal{L} = -\frac{1}{2} \sum_i \left( \frac{v_{\text{mod},i} - v_{\text{obs},i}}{\sigma_i} \right)^2 \quad (11)$$

where the summation $i$ is done over all the data points, $v_{\text{mod}}$ is the theoretical circular velocity from the MW models, $v_{\text{obs}}$ is the measured circular velocity, and $\sigma$ is the statistical uncertainty of the measurement (see Sect. 2.2).

<footnote n="3">³ We note that $R_\odot = 8.34$ kpc for this RC and 8.178 kpc for RC of Ou et al. (2023).</footnote>

<footnote n="4">⁴ In the present study, we accounted for the systematic uncertainties of Ou et al. (2023, see their Fig. 5) when deriving parameters from the corresponding RC. The values that cannot be seen in their Fig. 5 have been chosen to be 0.14.</footnote>

<footnote n="5">⁵ https://github.com/dfm/emcee</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** MCMC tests of Einasto dark matter parameters with several baryonic models (represented by coloured lines) for our RC (right panel) and that of Ou et al. (2023) RC.</figcaption></figure>

In the present study, we extrapolated the DM halo to the virial radius, denoted as $R_{vir}$, which encloses the virial or the DM contribution to the dynamical mass, $M_{vir}$. The virial radius is defined as the radius of a sphere within which the average density of dark matter is equal to 200 times the critical density of the Universe $\rho_{cr}$. Here, we adopted a critical density of $\rho_{cr} = 1.34 \times 10^{-7} M_{\odot} \text{pc}^{-3}$ (Hinshaw et al. 2013).

In order to properly estimate the dynamical mass, we investigated various baryonic models, as described in Sect. 3.1. The posterior distributions for both RC fits with different baryonic models are given in Fig. 6. We note that we converted the parameters of the virial DM mass $M_{vir}$, scale length $h$, and Einasto index $n$.

We calculated the dynamical mass as the sum of the baryonic mass and the DM halo virial mass. Results are presented in Table 4. We find a MW dynamical mass of $1.99_{-0.06}^{+0.09} \times 10^{11} M_{\odot}$ at $R_{vir} = 121.03_{-1.23}^{+1.80}$ kpc from the present work, which can be compared to the dynamical mass of $2.13_{-0.12}^{+0.17} \times 10^{11} M_{\odot}$ at $123.80_{-2.37}^{+3.21}$ kpc using the RC of Ou et al. (2023). The local DM density is found in the range of 0.011 to $0.012 M_{\odot} \text{pc}^{-3}$ (0.418–0.456 GeV cm$^{-3}$) for both RCs with different baryonic models. One may question the significance of such small error bars, given the fact that they lead to only very small values for the MW dynamical mass. On one hand, this is assuming that each point (and associated error bars) of the RC has been determined independently from the others, which is the underlying assumption of the MCMC method. This is likely true for the RC points, but not necessarily for the error bars, which may have led to an underestimation of the error bars on dynamical mass. On the other hand, error bars account for systematic uncertainties, and the latter may have been overestimated, as shown in Sect. 2.2.5.

We performed another test to evaluate the error bars of the MW dynamical mass. To this end, we arbitrarily replaced all the points of our RC by values provided by their upper error bars. By fitting these points with different baryonic models, we obtain dynamical masses ranging from $2.44$ to $2.53 \times 10^{11} M_{\odot}$.

The latter value would correspond to an absolute upper limit for the MW dynamical mass. However, Fig. 7 shows that Ou et al. (2023) RC points have larger error bars than in our study. Performing the same exercise with these values, one would find total mass values from $3.8$ to $5.4 \times 10^{11} M_{\odot}$. We consider the latter value as an absolute upper limit on the MW dynamical mass.

# 5. Discussion

## 5.1. The impact of asymmetric drift on the rotation curve measurement

We calculated the asymmetric drift following Eqs. (4.225), (4.227), and (4.228) (Strömberg's equation) in Binney & Tremaine (2008), and it could impact RCs. D23 provided velocities for RGB and OB stars that were not corrected for the asymmetric drift. Figure 7 compares the corresponding RCs to ours and that of Ou et al. (2023). This comparison shows that the non-corrected RGBs (magenta triangles) from D23 lag the corrected curves (see blue and red dots), which clearly illustrates the necessity of the asymmetric drift correction. OB stars (green triangles) are rotating faster, and as they rotate closer to the circular velocity, this would imply that both the correction in our study and that in Ou et al. (2023) are sufficient, at least for radii of greater than 10 kpc.

However, our correction does not seem entirely sufficient for $R < 10$ kpc, perhaps because of the unknown shape of the true density profile of the MW. This might also mean that the decrease in this region is even steeper than shown by the red and blue points. In the future, extending the very young star rotation curve towards inner disc regions, as well as to larger radii, may provide evidence supporting the decrease in the MW RC, and in particular might help to better determine the radius at which the decrease occurs.

Concerning systematic uncertainties introduced by the asymmetric drift, Ou et al. (2023, see their Fig. 5) showed that beyond 22 kpc, the asymmetric drift can contribute to relative systematic uncertainties at the level of 15%, which might cause the
<table><caption>**Table 4.** Estimated MW dynamical mass and associated Einasto profile parameters for the RC of the present work (TW) and for that of Ou et al. (2023, O23).</caption><thead><tr><th rowspan="2">Baryon model</th><th rowspan="2"><i>M</i><sub>bar</sub><br>[10<sup>11</sup> <i>M</i><sub>☉</sub>]</th><th colspan="2"><i>M</i><sub>dyn</sub><br>[10<sup>11</sup> <i>M</i><sub>☉</sub>]</th><th colspan="2"><i>M</i><sub>0</sub><br>[10<sup>11</sup> <i>M</i><sub>☉</sub>]</th><th colspan="2"><i>h</i><br>[kpc]</th><th colspan="2"><i>n</i></th></tr><tr><th>TW</th><th>O23</th><th>TW</th><th>O23</th><th>TW</th><th>O23</th><th>TW</th><th>O23</th></tr></thead><tbody><tr><td>B2</td><td>0.616</td><td>2.05<sup>+0.08</sup><sub>-0.06</sub></td><td>2.19<sup>+0.17</sup><sub>-0.12</sub></td><td>3.72<sup>+0.45</sup><sub>-0.70</sub></td><td>1.23<sup>+0.63</sup><sub>-0.58</sub></td><td>11.41<sup>+1.15</sup><sub>-1.62</sub></td><td>5.5<sup>+1.46</sup><sub>-1.56</sub></td><td>0.43<sup>+0.12</sup><sub>-0.09</sub></td><td>0.87<sup>+0.20</sup><sub>-0.15</sub></td></tr><tr><td>E dJ</td><td>0.607</td><td>1.97<sup>+0.09</sup><sub>-0.06</sub></td><td>2.07<sup>+0.15</sup><sub>-0.11</sub></td><td>3.72<sup>+0.36</sup><sub>-0.63</sub></td><td>1.82<sup>+0.64</sup><sub>-0.72</sub></td><td>12.30<sup>+1.10</sup><sub>-1.63</sub></td><td>7.3<sup>+1.46</sup><sub>-1.74</sub></td><td>0.40<sup>+0.13</sup><sub>-0.09</sub></td><td>0.73<sup>+0.18</sup><sub>-0.14</sub></td></tr><tr><td>E J</td><td>0.603</td><td>1.97<sup>+0.09</sup><sub>-0.06</sub></td><td>2.08<sup>+0.16</sup><sub>-0.11</sub></td><td>3.72<sup>+0.36</sup><sub>-0.70</sub></td><td>1.70<sup>+0.65</sup><sub>-0.72</sub></td><td>12.21<sup>+1.12</sup><sub>-1.68</sub></td><td>7.03<sup>+1.49</sup><sub>-1.79</sub></td><td>0.40<sup>+0.13</sup><sub>-0.09</sub></td><td>0.76<sup>+0.19</sup><sub>-0.14</sub></td></tr><tr><td>E CM</td><td>0.589</td><td>1.97<sup>+0.09</sup><sub>-0.06</sub></td><td>2.10<sup>+0.17</sup><sub>-0.12</sub></td><td>3.55<sup>+0.43</sup><sub>-0.73</sub></td><td>1.26<sup>+0.69</sup><sub>-0.61</sub></td><td>11.63<sup>+1.20</sup><sub>-1.77</sub></td><td>5.81<sup>+1.58</sup><sub>-1.71</sub></td><td>0.43<sup>+0.14</sup><sub>-0.10</sub></td><td>0.85<sup>+0.21</sup><sub>-0.16</sub></td></tr><tr><td>G dJ</td><td>0.575</td><td>1.98<sup>+0.09</sup><sub>-0.07</sub></td><td>2.14<sup>+0.17</sup><sub>-0.12</sub></td><td>3.47<sup>+0.51</sup><sub>-0.78</sub></td><td>1.02<sup>+0.64</sup><sub>-0.52</sub></td><td>11.34<sup>+1.28</sup><sub>-1.86</sub></td><td>4.99<sup>+1.52</sup><sub>-1.54</sub></td><td>0.45<sup>+0.15</sup><sub>-0.10</sub></td><td>0.94<sup>+0.21</sup><sub>-0.17</sub></td></tr><tr><td>G J</td><td>0.571</td><td>1.98<sup>+0.09</sup><sub>-0.06</sub></td><td>2.15<sup>+0.19</sup><sub>-0.13</sub></td><td>3.39<sup>+0.50</sup><sub>-0.76</sub></td><td>0.89<sup>+0.62</sup><sub>-0.51</sub></td><td>11.29<sup>+1.26</sup><sub>-1.82</sub></td><td>4.72<sup>+1.60</sup><sub>-1.64</sub></td><td>0.45<sup>+0.14</sup><sub>-0.10</sub></td><td>0.97<sup>+0.24</sup><sub>-0.18</sub></td></tr><tr><td>G CM</td><td>0.557</td><td>1.99<sup>+0.10</sup><sub>-0.07</sub></td><td>2.16<sup>+0.18</sup><sub>-0.13</sub></td><td>3.31<sup>+0.58</sup><sub>-0.91</sub></td><td>0.64<sup>+0.56</sup><sub>-0.38</sub></td><td>10.63<sup>+1.40</sup><sub>-2.01</sub></td><td>3.85<sup>+1.52</sup><sub>-1.35</sub></td><td>0.49<sup>+0.16</sup><sub>-0.11</sub></td><td>1.06<sup>+0.23</sup><sub>-0.19</sub></td></tr></tbody><tfoot><tr><td>Average</td><td></td><td>1.99<sup>+0.09</sup><sub>-0.06</sub></td><td>2.13<sup>+0.17</sup><sub>-0.12</sub></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tfoot></table>

peculiar point around 23 kpc. Figure 8 of Wang et al. (2023a) shows that the declining slope is shared for stars selected at different heights, for which differences are most likely caused by the asymmetric drift effect. Points near and above $R = 24$ kpc might also be affected by the asymmetric drift systematic errors, which is similar to the findings of Ou et al. (2023) for their RC.

## 5.2. Does the Milky Way warp have a significant impact on the rotation curve measurement?

The MW has a warped disc. Observations show that the MW disc is flat out to roughly the Solar Circle, where it then bends upwards and downwards in the northern and southern hemispheres. The amplitude of the warp clearly increases with radius and varies with azimuthal angle. The position and kinematics of the Galactic warp also depend on the stellar populations (Poggio et al. 2018; Wang et al. 2020; Chrobáková et al. 2020a; Li et al. 2023).

Wang et al. (2020) analysed the stellar warp of different stellar populations by combining LAMOST DR4 and Gaia DR2. These authors concluded that the Galactic warp might not be caused by the gravitational interaction scenarios but by the gas infall process, and gravitational interactions such as that due to the Sgr passages are also addressed.

On the other hand, there has been no major merger in the MW for between 9 and 10 Gyr, which corresponds to the GSE event, as described in Sect. I. Interactions with nearby dwarfs or satellites are not expected to have a great impact on disc stability, though Sgr passages near the disc edge have been suggested to have affected or even possibly formed the warp (Bailin 2003). Even in the latter case, if the Sgr main interaction with the disc occurred more than 5 Gyr ago, it is likely that the MW disc is in relative equilibrium (quasi-equilibrium within 30 kpc)<fnref n="6" /> with some oscillations. Those oscillating asymmetries on either side of the galactic plane are called corrugations. From a sample of 40 nearby low-inclination disc galaxies, Urrejola-Mora et al. (2022) identified that 20% of the galaxies exhibit vertically perturbed galactic discs, which could be described by corrugations. One of the most famous corrugations in the MW is

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Circular velocity of the Milky Way. The red (blue) data points are the measurements computed in the present work (and in Ou et al. 2023), with error bars including systematic uncertainties. These are compared with data points from D23 (OB: magenta triangles, and RGB: green triangles, stars) that have not been corrected for asymmetric drift.</figcaption></figure>

the Monoceros Ring located at low Galactic latitudes. Although this is more likely to be explained by disc flaring (Wang et al. 2018; Bergemann et al. 2018), scenarios of a perturbed disc by an ancient disrupted dwarf galaxy have also been proposed (Conn et al. 2007, 2012; Johnston et al. 2017). Corrugations start from the Sun and at least four ripples are seen in the disc outskirts (Newberg et al. 2002; Ibata et al. 2003; Xu et al. 2015); these blur our representation of a flat disc and can produce wiggles in the RC that are different from those due to spiral arms.

Chrobáková et al. (2020b) used N-body simulations and found that the Jeans equation could provide a reasonable approximation to the system dynamics if the amplitude of the radial velocity component is significantly smaller than the azimuthal one. It is possible to assume an axisymmetric steady disc to measure the RC globally using the Jeans equation. In Fig. 1 (see Appendix A), we also compared the RC with different limits for the vertical heights and they show good consistency even in the outer disc. The warp could affect spatial and velocity measurements, and especially the vertical velocities, which are expected to be its main kinematic signature. In order to measure the effect of the vertical velocity, we analysed the neglected cross-term

<footnote n="6">⁶ At 30 kpc, a star rotating at 170 km s⁻¹ would have the time to make 4.5 orbits, i.e. a sufficiently large number to consider the system to be at equilibrium (Gnedin & Ostriker 1999).</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 8. Enclosed dynamical mass of the MW as a function of radius.</figcaption></figure>

$\langle V_R V_z \rangle$ in detail (see Sect. 2.2.1). We also note that Wang et al. (2023a) applied LIM to the full Gaia DR3 data set with different stellar populations, which helps to reduce the warp impact because we know that the warp is much stronger for young populations.

### 5.3. Have we reached the Keplerian decline of the Milky Way rotation curve?

Our RC measurements based on Gaia DR3 (as well as that from Ou et al. 2023) demonstrate its significant decline. Additionally, the slope of the declining RC is steeper compared to the previous studies (Eilers et al. 2019; Mróz et al. 2019) based on Gaia DR2. This also implies that the increase in cumulative mass at larger radii is minimal. We apply an approximate conversion from the measured circular velocity to the enclosed dynamical mass using the following relation:

$$M_{dyn,i} = \frac{v_{circ,i}^2 R_i}{G}, \qquad (12)$$

where G is the gravitational constant and $v_{circ,i}$ is the circular velocity at each radius $R_i$. The conversion is not accurate in a non-spherical potential but the difference is very small at large radii, where the spherical DM component dominates (see Fig. 5). By applying Eq. (12), we determined that beyond R > 19 kpc, the enclosed mass barely varies (see Fig. 8), remaining within 1.9 to 2.0 × 10¹¹ M$_{\odot}$, which is remarkably close to our estimate made at a significantly larger radius (see Table 4). The small decay of enclosed dynamical mass at large radii (R > 23 kpc) in Fig. 8 cannot be physical. However, given that the amplitude of the decay is much smaller than the error bars, it has no incidence on the validity of this work.

Given that the RC drops faster at large radii and because of the agreement between enclosed and virial masses, we conducted a test to assess whether or not the MW RC has reached a Keplerian decline beyond R > 19 kpc. In Fig. 9, we compare the Keplerian decline using the same enclosed mass at each radius beyond R > 19 kpc, which helps us to verify whether or not the MW RC at the outskirts can be fitted this way. Figure 9 displays the best fit for a Keplerian decline to the RCs from our measurement (top panel) and that of Ou et al. (2023, bottom panel), respectively. At large radii, this fits well with the two estimated RCs, suggesting that stellar rotation at the outskirts follows a Keplerian decline, which implies that the enclosed mass within a certain radius is sufficient to account for the observed velocities, without requiring an increase in mass at

<figure><img src="image_10.png" /><figcaption>Fig. 9. Rotation curve with the best fit of the Keplerian decline for our RC measurement (top panel) and that from Ou et al. (2023, bottom panel).</figcaption></figure>

larger radii. Interestingly, we find a consistent enclosed mass for the Keplerian decline in both RC measurements, which amounts to approximately 1.95 × 10¹¹ M$_{\odot}$.

In order to further test for the presence of a Keplerian decline at large radii, we assumed a circular velocity profile following:

$$V(R) = AR^{\gamma}, \qquad (13)$$

where A is an amplitude parameter, and γ is the exponential slope of the outer RC. The posterior distributions are presented in Fig. 10. By sampling γ using a flat prior between -10 and 5, we found that the slope γ beyond R > 19 kpc is -0.47+0.15 and -0.56+0.23 for the RC of this work and that of Ou et al. (2023), respectively. This suggests that the MW RC is consistent with a Keplerian decline (γ = -0.5) at large radii, while a flat RC is rejected at a 3σ significance level (see green arrows in Fig. 10).

### 5.4. Comparison of the MW dynamical mass with that based on different mass tracers

The RC of the Milky Way is the most accurate tracer for estimating the enclosed dynamical mass within a range of approximately 30 kpc, because rotating stars in the disc are likely at equilibrium with the total potential associated with the enclosed dynamical mass. We note that our dynamical mass estimate is primarily constrained by the RC within the range of R = 9–27 kpc and is then extrapolated to the virial radius. The extrapolation is unlikely to bring an additional mass component, because the MW RC seems to follow a Keplerian decline at large radii.

Other estimates derived from stellar streams (Vasiliev et al. 2021; Koposov et al. 2023), globular clusters (Eadie & Jurić 2019; Posti & Helmi 2019), the Magellanic Cloud (Correa Magnus & Vasiliev 2022), or dwarf galaxies (Cautun et al. 2020; Li et al. 2020; Slizewski et al. 2022) provide virial masses
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Posterior distributions for the parameters *A* and *γ*. We show fits to the RC of this work (left panel) and to that of Ou et al. (2023, right panel) including systematics uncertainties. Histograms of parameters *A* and *γ* correspond to the maximum density distribution. Red solid lines represent the median of the posterior distributions. Black dashed lines correspond to 1σ uncertainty. Blue dashed lines correspond to the slope of the Keplerian decline, *n* = −0.5. Green arrows indicate the slope of a flat RC (*n* = 0).</figcaption></figure>

ranging from 7 to $11 \times 10^{11} M_{\odot}$ at about 150 to 250 kpc. Our predicted virial mass from the MW RC is then considerably lower than these previous estimates.

However, the findings of recent studies of the orbital histories of halo inhabitants may challenge our previous understanding. Hammer et al. (2023) found that most dwarf galaxies are relative newcomers. This is due to the expected linear relationship between infall times and the logarithm of the total orbital energy (Rocha et al. 2012). As dwarf-galaxy orbital energies are larger than those of former events such as the GSE and the Sgr infall, they likely entered the MW halo less than 3 Gyr ago. Moreover, Li et al. (2021) found that the dwarf galaxies are highly concentrated near their pericentre. Consequently, the assumption of virial equilibrium for dwarf-galaxy orbital motions is likely invalid, and it is not surprising that one finds very large mass values when using these orbital motions as mass tracers for the MW. This is well illustrated by Boylan-Kolchin et al. (2013) who showed that assuming Leo I as a bound satellite would lead to a significant overestimate of the MW mass.

After using both globular clusters and the MW RC from Gaia DR2, Wang et al. (2022b) found that excluding Crater and Pyxis, which possess large amounts of orbital energy, leads to a decrease in the MW mass estimate from $5.73_{-0.58}^{+0.76} \times 10^{11} M_{\odot}$ to $5.36_{-0.68}^{+0.81} \times 10^{11} M_{\odot}$ when assuming an Einasto profile.

These recent findings highlight the fact that our understanding of the MW structure is changing. They also emphasise the need to verify whether or not various tracers have their orbital velocity in equilibrium with the MW potential in order to avoid systematic overestimates of the total mass.

### 5.5. Is our Galaxy exceptional or is its rotation curve related to the Gaia methodology used to recover it?

We investigated the literature on searching for declining spiral galaxy RCs to find out whether or not some of them

may show a Keplerian decline. In reviewing galaxy discs, van der Kruit & Freeman (2011) mentioned that none spiral galaxies show a decline in their RCs, which can be associated with a cut-off in the mass distribution, so that in no case has the RC been traced to the limit of the dark matter distribution. In their analysis of well-studied spirals from the THINGS project, de Blok et al. (2008) found that a declining RC such as that of M81 is likely caused by galaxy interactions.

Noordermeer et al. (2007, see also Zobnina & Zasov 2020) specifically studied spirals with a declining RC and found no RCs with a fully Keplerian decline in the outer regions, indicating that we have not yet reached the point where the mass density becomes negligible, except perhaps for UGC 4458 (NGC 2599). However, the latter shows a bright UV nucleus that may be active, and its disc morphology also appears quite different from that of the MW. Dicaire et al. (2008) identified NGC 7793 as a possible candidate for exhibiting a Keplerian decline; its DM content is surprisingly low for a dwarf galaxy, and is smaller than that of stars at all radii.

The MW may indeed be relatively exceptional if it is the only isolated grand-design spiral showing a Keplerian decline in its RC, which is consistent with a possible cut-off of its mass distribution. One might wonder whether or not this declining RC is related to its four-arm structure discovered by Georgelin & Georgelin (1976). Four-armed spirals are rare in nearby Milky Way-type galaxies, although galaxies with two arms frequently exhibit bifurcations resulting in several arm segments, which could also be the case for the Milky Way.

Alternatively, our Galaxy may be exceptional due to its relatively quiet past history (Hammer et al. 2007; Belokurov et al. 2018; Haywood et al. 2018), having experienced no major merger for ~9 Gyr. We must also consider the possibility that the methodology used by Gaia to recover its full 6D space-velocity parameters for very large numbers of stars may be a contributing factor to the exceptional nature of the RC we obtain based on
these data; it certainly contrasts with the less constrained RC of external galaxies.

## 6. Conclusions

Gaia DR3 has led to significant progress in our capacity to estimate the MW RC when compared to Gaia DR2. Three different studies have derived the RC of the MW: One used the whole and very large Gaia DR3 catalogue with distances estimated from parallaxes, and then averaged into specific 6D cells using LIM (Wang et al. 2023a). The second (Ou et al. 2023) is based on a smaller number of RGB stars (33 335), for which distances are estimated from spectrophotometry. The third (Zhou et al. 2023) is based on 58 000 bright RGB stars, and in the present study we re-evaluated their distance estimates, rendering their RC consistent with that of the two other studies.

Here, we carried out a full analysis of the systematic uncertainties that can affect the Wang et al. (2023a) study. We then compared the different RCs from Gaia DR3, and find that they have consistent rotational velocity values from 9 to 27 kpc. This indicates a robust and significant decline of the MW RC, with a slope of $\beta = -(2.18 \pm 0.23) \text{ km s}^{-1} \text{ kpc}^{-1}$, which is based on combining the very similar slopes from this study and that from Ou et al. (2023). The decrease in velocity between 19.5 and 26.5 kpc is approximately 30 km s$^{-1}$.

The estimated MW dynamical mass is consistent with $1.99_{-0.06}^{+0.09} \times 10^{11} M_{\odot}$ at 121.03$^{+1.80}_{-1.23}$ kpc (from the RC of this study) and with $2.13_{-0.12}^{+0.17} \times 10^{11} M_{\odot}$ at 123.80$^{+3.21}_{-2.37}$ kpc (from the RC of Ou et al. 2023, respectively)$^7$. The relatively small size of the error bars is perhaps a result of the assumption that the two RC points and their associated error bars were independently determined.

Consequently, the ratio of DM to baryonic mass is only a factor of about 3, instead of a factor of approximately 6 from $\Lambda$CDM (Planck Collaboration VI 2020), which suggests that baryons are not missing in our Galaxy. A small dynamical mass for the MW may also impact mass estimations for the LMC, which is mostly constrained by its ratio to that of the MW; for example, when studying the induced sloshing of the Galactic halo (Erkal et al. 2021; Conroy et al. 2021). If the total MW mass is as small as $2.06 \times 10^{11} M_{\odot}$, the LMC total mass would be from 2 to $3 \times 10^{10} M_{\odot}$. Interestingly, the latter value is consistent with the modelling of the Magellanic Stream through ram pressure as shown by Hammer et al. (2015) and Wang et al. (2019, see also LMC mass predictions from Wang et al. 2022a), which would resolve the numerous difficulties in reproducing it with tidal tail models (Besla et al. 2012; Lucchini et al. 2020, 2021).

We conclude that the MW RC cannot be consistent with a flat RC at a significance of $3\sigma$, and that our findings suggest a Keplerian decline occurring at radii of greater than 19 kpc. The Keplerian decline indicates the point where the mass density becomes negligible (Zobnina & Zasov 2020). Some spiral galaxies present a declining RC, but at large radii they appear to flatten out, meaning that their RC has not been traced to the outermost extent of the dark matter distribution (de Blok et al. 2008; van der Kruit & Freeman 2011). This contrasts with the MW, whose accretion history (Hammer et al. 2023, and references therein) shows no major merger for 8 to 10 Gyr, while half of the spiral galaxies underwent their last major merger more recently (Hammer et al. 2005, 2009; Puech et al. 2012). It would

be interesting to study the impact of relatively recent assembly events on the RC at the outskirts of spiral galaxies.

In many respects, a Keplerian decline for the MW RC may appear quite exceptional. This could be due to the extraordinarily quiet history of our Galaxy, or to the very different methodology used by Gaia to calculate its kinematics compared to that used to study external galaxies.

**Acknowledgements.** We thank the anonymous referee for their comments and suggestions. Y.-J.J. acknowledges financial support from the China Scholarship Council (CSC) No. 202108070090. J.-L.W. is also grateful to be supported by the CSC No. 202210740004. L.C. acknowledges the Chilean Agencia Nacional de Investigación y Desarrollo through the grant Fondecyt Regular 1210992. We are grateful to Chrobáková Žofia, López-Corredoira Martín and Francesco Sylos Labini for their invaluable assistance in producing the datasets and tests of Wang et al. (2023a). We thank Yang Huang for sharing their data from Zhou et al. (2023). We are grateful for the support of the International Research Program Tianguan, which is an agreement between the CNRS in France, NAOC IHEP, and the Yunnan Univ. in China. The data underlying this article will be shared on request to the corresponding author.

## References

Babcock, H. W. 1939, Lick Obs. Bull., 498, 41

Bailin, J. 2003, ApJ, 583, L79

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Bergemann, M., Sesar, B., Cohen, J. G., et al. 2018, Nature, 555, 334

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2012, MNRAS, 421, 2109

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press)

Bosma, A. 1978, PhD Thesis, University of Groningen, The Netherlands

Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115

Boylan-Kolchin, M., Bullock, J. S., Sohn, S. T., Besla, G., & van der Marel, R. P. 2013, ApJ, 768, 140

Calchi Novati, S., & Mancini, L. 2011, MNRAS, 416, 1292

Cautun, M., Benítez-Llambay, A., Deason, A. J., et al. 2020, MNRAS, 494, 4291

Chemin, L., de Blok, W. J. G., & Mamon, G. A. 2011, AJ, 142, 109

Chrobáková, Ž., Nagy, R., & López-Corredoira, M. 2020a, A&A, 637, A96

Chrobáková, Ž., López-Corredoira, M., Sylos Labini, F., Wang, H. F., & Nagy, R. 2020b, A&A, 642, A95

Conn, B. C., Lane, R. R., Lewis, G. F., et al. 2007, MNRAS, 376, 939

Conn, B. C., Noël, N. E. D., Rix, H.-W., et al. 2012, ApJ, 754, 101

Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Nature, 592, 534

Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610

de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648

de Jong, J. T. A., Yanny, B., Rix, H.-W., et al. 2010, ApJ, 714, 663

de Salas, P. F., Malhan, K., Freese, K., Hattori, K., & Valluri, M. 2019, JCAP, 2019, 037

Dicaire, I., Carignan, C., Amram, P., et al. 2008, AJ, 135, 2038

Eadie, G., & Jurić, M. 2019, ApJ, 875, 159

Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, ApJ, 871, 120

Einasto, J. 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Erkal, D., Deason, A. J., Belokurov, V., et al. 2021, MNRAS, 506, 2677

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Gaia Collaboration (Drimmel, R., et al.) 2023a, A&A, 674, A37

Gaia Collaboration (Vallenari, A., et al.) 2023b, A&A, 674, A1

Georgelin, Y. M., & Georgelin, Y. P. 1976, A&A, 49, 57

Gnedin, O. Y., & Ostriker, J. P. 1999, ApJ, 513, 626

Goodwin, S. P., Gribbin, J., & Hendry, M. A. 1998, The Observatory, 118, 201

Hammer, F., Flores, H., Elbaz, D., et al. 2005, A&A, 430, 115

Hammer, F., Puech, M., Chemin, L., Flores, H., & Lehnert, M. D. 2007, ApJ, 662, 322

Hammer, F., Flores, H., Puech, M., et al. 2009, A&A, 507, 1313

Hammer, F., Yang, Y. B., Flores, H., Puech, M., & Fouquet, S. 2015, ApJ, 813, 110

Hammer, F., Li, H., Mamon, G. A., et al. 2023, MNRAS, 519, 5059

Haywood, M., Di Matteo, P., Lehnert, M. D., et al. 2018, ApJ, 863, 113

Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85

Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19

Hogg, D. W., Eilers, A.-C., & Rix, H.-W. 2019, AJ, 158, 147

Ibata, R. A., Irwin, M. J., Lewis, G. F., Ferguson, A. M. N., & Tanvir, N. 2003, MNRAS, 340, L21

<footnote n="7">⁷ We note that the two estimates are very similar and the average value would be $2.06_{-0.13}^{+0.24} \times 10^{11} M_{\odot}$.</footnote>
Iocco, F., Pato, M., & Bertone, G. 2015, Nat. Phys., 11, 245

Jiao, Y., Hammer, F., Wang, J. L., & Yang, Y. B. 2021, A&A, 654, A25

Johnston, K. V., Price-Whelan, A. M., Bergemann, M., et al. 2017, Galaxies, 5, 44

Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864

Karukes, E. V., Benito, M., Iocco, F., Trotta, R., & Geringer-Sameth, A. 2020, JCAP, 2020, 033

Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260

Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936

Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80

Laporte, C. F. P., Johnston, K. V., Gómez, F. A., Garavito-Camargo, N., & Besla, G. 2018, MNRAS, 481, 286

Li, Z.-Z., Qian, Y.-Z., Han, J., et al. 2020, ApJ, 894, 10

Li, H., Hammer, F., Babusiaux, C., et al. 2021, ApJ, 916, 8

Li, X., Wang, H.-F., Luo, Y.-P., et al. 2023, ApJ, 943, 88

Lindegren, L., Bastian, U., Biermann, M., et al. 2021a, A&A, 649, A4

Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021b, A&A, 649, A2

Lucchini, S., D'Onghia, E., Fox, A. J., et al. 2020, Nature, 585, 203

Lucchini, S., D'Onghia, E., & Fox, A. J. 2021, ApJ, 921, L36

Lucy, L. B. 1974, AJ, 79, 745

Lundmark, K. 1925, MNRAS, 85, 865

Mayall, N. U. 1951, Pub. Mich. Obs., 10, 19

Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245

Noordermeer, E., van der Hulst, J. M., Sancisi, R., Swaters, R. S., & van Albada, T. S. 2007, MNRAS, 376, 1513

Ou, X., Eilers, A. C., Necib, L., & Frebel, A. 2023, ArXiv e-prints [arXiv:2303.12838]

Planck Collaboration VI. 2020, A&A, 641, A6

Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21

Posti, L., & Helmi, A. 2019, A&A, 621, A56

Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66

Puech, M., Hammer, F., Hopkins, P. F., et al. 2012, ApJ, 753, 128

Queiroz, A. B. A., Anders, F., Chiappini, C., et al. 2023, A&A, 673, A155

Retana-Montenegro, E., van Hese, E., Gentile, G., Baes, M., & Frutos-Alfaro, F. 2012, A&A, 540, A70

Rocha, M., Peter, A. H. G., & Bullock, J. 2012, MNRAS, 425, 231

Rubin, V. C., Ford, W. K., Jr., & Thonnard, N. 1978, ApJ, 225, L107

Sackett, P. D. 1997, ApJ, 483, 103

Schmidt, M. 1965, in Galactic Structure, eds. A. Blaauw, & M. Schmidt (Chicago: University of Chicago Press), 513

Slizewski, A., Dufresne, X., Murdock, K., et al. 2022, ApJ, 924, 131

Sofue, Y., Honma, M., & Omodaka, T. 2009, PASJ, 61, 227

Stanek, K. Z., Udalski, A., Szymański, M., et al. 1997, ApJ, 477, 163

Sylos Labini, F., Chrobáková, Ž., Capuzzo-Dolcetta, R., & López-Corredoira, M. 2023, ApJ, 945, 3

Urrejola-Mora, C., Gómez, F. A., Torres-Flores, S., et al. 2022, ApJ, 935, 20

van der Kruit, P. C., & Freeman, K. C. 2011, ARA&A, 49, 301

Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279

Wang, J., Shi, J., Pan, K., et al. 2016, MNRAS, 460, 3179

Wang, H.-F., Liu, C., Xu, Y., Wan, J.-C., & Deng, L. 2018, MNRAS, 478, 3367

Wang, J., Hammer, F., Yang, Y., et al. 2019, MNRAS, 486, 5907

Wang, H. F., López-Corredoira, M., Huang, Y., et al. 2020, ApJ, 897, 119

Wang, J., Hammer, F., & Yang, Y. 2022a, MNRAS, 515, 940

Wang, J., Hammer, F., & Yang, Y. 2022b, MNRAS, 510, 2242

Wang, H.-F., Chrobáková, Ž., López-Corredoira, M., & Sylos Labini, F. 2023a, ApJ, 942, 12

Wang, J., Cao, Z., Huang, Y., & Yuan, H. 2023b, Res. Astron. Astrophys., 23, 025020

Wegg, C., Gerhard, O., & Portail, M. 2016, MNRAS, 463, 557

Xu, Y., Newberg, H. J., Carlin, J. L., et al. 2015, ApJ, 801, 105

Zhou, Y., Li, X., Huang, Y., & Zhang, H. 2023, ApJ, 946, 73

Zobnina, D. I., & Zasov, A. V. 2020, Astron. Rep., 64, 295
**Appendix A: The impact of the z selection on the**
**RC between 9 and 13 kpc**

Figure A.1 indicates that by limiting the Wang et al. (2023a) data to |z| < 1 kpc, their RC from 9 to 13 kpc becomes consistent with that of Ou et al. (2023).

<figure><img src="image_12.png" /><figcaption>**Fig. A.1.** Comparison of the RC using different selections in the vertical direction.</figcaption></figure>
**Appendix B: Comparison of distance estimates**

Figure B.1 compares the distance estimates by Zhou et al. (2023) with those of three different methods on the APOGEE data. One is a modified version of Wang et al. (2016), in which Gaia EDR3 parallax has been incorporated as a prior to constrain stellar distances as described in Wang et al. (2023b). Other distance estimates are from Hogg et al. (2019), which were used by Eilers et al. 2019), and from the results of StarHorse

(Queiroz et al. 2023). This comparison shows that Zhou et al. (2023) always overestimate distances for $R > 10$ kpc stars. To evaluate the consequences of such a bias, Figure B.2 provides a very rough analysis by correcting each Zhou et al. (2023) RC point accordingly, which suffices to reconcile their RC with those from this study and from Ou et al. (2023). However, a better analysis, where the distances are corrected for individual stars, is still needed.

<figure><img src="image_13.png" /></figure>

**Fig. B.1.** Comparison of distance estimates made by Zhou et al. (2023) to those made by Wang et al. (2016), Hogg et al. (2019), and StarHorse (Queiroz et al. 2023) in the top panels, respectively. The solid line presents a one-to-one correspondence. The dashed line presents the best linear fit on the logarithmic scale. Red points and error bars indicate the median and standard deviation per 400 data points. The ratios $\Delta d/d_{zhou}$ are shown in the bottom panels.

<figure><img src="image_14.png" /></figure>

**Fig. B.2.** Comparison of RCs after a rough correction of the distance applied to the points from Zhou et al. (2023, see text).