Mapping the Milky Way Disk with Gaia DR3: 3D Extended Kinematic Maps and Rotation Curve to ≈30 kpc

The Gaia mission (Gaia Collaboration et al. 2016) is providing the most detailed Milky Way survey to date by measuring stellar astrometry, photometry, and spectroscopy. From this data set it is possible to derive the spatial distribution, kinematics, and many other physical properties of the Milky Way. The third Gaia data release (DR3; Gaia Collaboration et al. 2022b) contains the same photometric magnitude and astrometric information as the previous Gaia data EDR3 releases but for a wider sample of stars with new determinations of spectra, radial velocity, chemical abundance, value-added catalogs, etc.

Concerning the kinematics of our Galaxy, the Gaia data sets provide the full astrometric solution—positions on the sky, parallax, and proper motion. Gaia DR3 has also provided a significant increase of the stars' line-of-sight velocity catalog, from 7,209,831 in DR2 to 33,812,183 in DR3 (Katz et al. 2022).

Since the first data release, Gaia has played a major part in revealing the kinematics of our Galaxy. Antoja et al. (2016, 2017) used Gaia DR1 and mock data to probe the influence of the spiral arms on Galactic kinematics and found that the typical difference in transverse velocity at symmetric longitudes is about 2 km s⁻¹, but can be larger than 10 km s⁻¹ at some longitudes and distances.

By using the DR2, Gaia Collaboration et al. (2018) revealed streaming motions in all three velocity components and found that the vertical velocities show a superposition of modes. The same work found small-amplitude perturbations in the velocity dispersion as well.

Kawata et al. (2018) have carried out a kinematic analysis for radii R < 13 kpc, a range in which the relative error in the distance was lower than 20%. Poggio et al. (2018) extended the analyzed region up to R < 15 kpc in combination with Two Micron All Sky Survey (2MASS) photometry, with the aim of characterizing vertical motions. By using 5 million stars from the Gaia DR2 catalog belonging to the Milky Way disk, Ramos et al. (2018) have also detected many kinematic asymmetries whose azimuthal velocity decreased with the galactic radius.

López-Corredoira & Sylos Labini (2019) have presented extended kinematic maps of the Galaxy with Gaia DR2, including the region where the relative error in distance was between 20% and 100%. That was possible thanks to the use of a statistical deconvolution algorithm of the parallax errors named Lucy's inversion method (LIM; Lucy 1974 and references therein).

By applying this method to the Gaia DR2 data set and including line-of-sight velocity measurements, López-Corredoira & Sylos Labini (2019) have extended the range of distances for the kinematic analyses by ≈7 kpc with respect to those presented by Gaia Collaboration et al. (2018), thus adding the range of Galactocentric distances between 13 and 20 kpc to the previous maps. They found velocity gradients of about 40 km s⁻¹ in the radial and azimuthal directions and 10 km s⁻¹ in the vertical direction, as well as north–south asymmetries.

Bennett & Bovy (2019) made detailed measurements of the wave-like north–south asymmetry in the vertical stellar counts, which showed some deficits at heights ≈0.4, ≈0.9, and ≈1.5 kpc and peaks at ≈0.2, ≈0.7 and ≈1.1 kpc. The mean vertical velocity is also found to show a north–south symmetric dip at ≈0.5 kpc with an amplitude of ≈2 km s⁻¹. By using OB type stars samples and red giant branch (RGB) samples, Romero-Gómez et al. (2019) found that the median vertical proper-motion values show a clear vertical modulation toward the anticenter.

Antoja et al. (2021) used Gaia EDR3 to reveal that oscillations in median rotational and vertical velocities vary with radius, disk stars having large angular momentum moving vertically upward, and disk stars having slightly lower angular momentum moving preferentially downward.

Drimmel et al. (2022) traced the spiral structure and kinematics using the young stellar population from the Gaia DR3 data set and found that the local arms are at least 8 kpc long, the outer arms are consistent with those seen in the Hi survey, and the Perseus arms continue toward the third quadrant. In addition, by analyzing the RGB samples with measured velocities, they found the streaming motions in the outer disk that may be associated with the spiral arms or bar dynamics.

Our aim in this paper is to derive kinematic maps of the Galaxy up to R ≈ 30 kpc by applying LIM on the Gaia DR3 data set only from the observational points of view. There are two advantages passing from DR2 to DR3: on the one hand, the subsample including radial velocities extracted from Gaia DR3 has a fainter magnitude limit, i.e., G_RVS = 14 (Katz et al. 2022), than that from Gaia DR2 (G_RVS = 12) and thus many more sources; on the other hand, parallax errors become lower passing from DR2 to DR3. Indeed, according to Lindegren et al. (2021), uncertainties are smaller in DR3 with respect to DR2 by a factor 0.80 for parallaxes and positions and by a factor 0.51 for proper motions.

Moreover, by assuming that the Galaxy is in a steady state and that the gravitational potential is axisymmetric, i.e., by neglecting perturbations, we then can determine, by means of the Jeans equation, the velocity rotation curve V_C (R) or V_ϕ up to ≈30 kpc. This allows us to confirm the recent results by Eilers et al. (2019), who found that the V_C (R) is linearly declining up to 25 kpc with a slope of ≈−1.7 ± 0.1 km s⁻¹ kpc⁻¹ (systematic uncertainty of 0.46 km s⁻¹ kpc⁻¹). In addition, we present estimations of the rotation curve at different heights and compare them with the results by Chrobáková et al. (2020); most notably, we find, at small radii, a marked dependence of V_C (R) on Z, which will be shown below.

The paper is structured as follows: In Section 2 we briefly describe the Gaia DR3 sample we use in this work. Then, in Section 3 we outline essential elements of LIM, while in Section 4 we present results about the 3D kinematics and the rotation curve. Finally, in Section 5 we draw our main conclusions.

Gaia DR3 (Gaia Collaboration et al. 2022b) provides line-of-sight velocities of more than 33 million sources with a limiting magnitude of G ≈ 14. The full astrometric solution has been done as five-parameter fit for 585 million sources and a six-parameter fit for 882 million sources. The median parallax uncertainty is 0.01−0.02 mas for G < 15, 0.05 mas at G = 17, 0.4 mas at G = 20, and 1.0 mas at G = 21. The uncertainty in the determination of the proper motion is 0.02−0.03 mas yr⁻¹ for stars with G < 15, 0.07 mas yr⁻¹ for stars with G = 17, 0.5 mas yr⁻¹ for stars with G = 20, and 1.4 mas yr⁻¹ for stars with G = 21 (Lindegren et al. 2021).

Gaia DR3 also contains a release of magnitudes estimated from the integration of Radial Velocity Spectrometer (RVS) spectra for a sample of about 32.2 million stars brighter than G_RVS ≈ 14 (or G ≈ 15). Sartoretti et al. (2022) have described the data used and the approach adopted to derive and validate the G_RVS magnitudes published in DR3. They also provide estimates of the G_RVS passband and associated G_RVS zero-point. Recio-Blanco et al. (2022) have summarized the stellar parameterization of the Gaia RVS spectra (around 5 million) performed by the GSP-Spec module and published it as part of Gaia DR3. The formal median precision of radial velocities is 1.3 km s⁻¹ at G_RVS = 12 and 6.4 km s⁻¹ at G_RVS = 14. The velocities zero-point exhibit a small systematic trend with magnitude starting around G_RVS = 11 and reaching about 400 m s⁻¹ at G_RVS = 14; to take into account this trend, a correction formula is provided. Note that the Gaia DR3 velocity scale is in satisfactory agreement with APOGEE, GALAH, GES, and RAVE (Katz et al. 2022).

The Gaia DR3 data set provides for each star the parallax π, the Galactic coordinates (ℓ, b), the line-of-sight velocity V_r , and the two proper motions in equatorial coordinates ${\mu }_{\alpha }\cos \delta$ and μ_δ . These six variables allow us to place each star in the 6D phase space, i.e., 3D spatial coordinates plus 3D velocity coordinates. The Galactocentric position in cylindrical coordinates has components R, the Galactocentric distance, ϕ, the Galactocentric azimuth, and Z, the vertical distance. The Galactocentric velocity in cylindrical coordinates has components V_R , the radial velocity, V_ϕ , the azimuthal velocity, and V_Z , the vertical velocity. Details to transform respectively (π, ℓ, b) into Galactocentric Cartesian (X, Y, Z) or cylindrical (R, ϕ, z) coordinates and $({V}_{r},{\mu }_{\alpha }\cos \delta ,{\mu }_{\delta })$ in (V_R , V_ϕ , V_Z ) are given in López-Corredoira & Sylos Labini (2019). The 3D velocity we used is obtained by assuming that the location of the Sun is R_⊙ = 8.34 kpc (Reid et al. 2014) and Z_⊙ = 27 pc (Chen et al. 2001). We use the solar motion values [U_⊙, V_⊙, W_⊙] = [11.1, 12.24, 7.25] km s⁻¹ (Reid et al. 2014). The value of the circular speed of the LSR is 238 km s⁻¹ (Schönrich et al. 2010). Note that different solar values have negligible effect on this work.

It is known that Gaia parallax uncertainties increase with distance, in particular beyond 5 kpc from the Sun, and that they are also dependent on the star's magnitude. The large dispersion of parallaxes corresponds to a large dispersion of stellar distances: such dispersion affects the value of mean distance at large heliocentric distances, which is indeed overestimated in a direct measurement with respect to the real one. In order to solve the problem of the deconvolution of the Gaussian errors with large rms values and the asymmetric parallax uncertainties, we adopt LIM. For more details we refer the reader to López-Corredoira & Sylos Labini (2019), where Monte Carlo simulations to test such a technique are also presented.

The basic idea of the method is simple: in order to explore Galactic regions where relative parallax errors are larger than 20%, i.e., for R > 18 kpc, we can reduce the error in the distance determination by making averages over many stars. We thus divide the observed Galactic region into N_cells cells, each containing many stars. We want to determine the average value of the velocity components and their dispersion in each cell. LIM provides the deconvolved distribution of sources along the line of sight; for this quantity we can estimate the mean heliocentric distance of all the stars included in a bin centered at parallax π and the corresponding variance ${\sigma }_{\overline{r}}^{2}(\pi )$. It is worth noticing that LIM is model independent: we can recover information on the stellar distribution without introducing any prior. In addition, LIM does not provide the distance for each star; rather, we only have statistical determination of the distance of a 3D cell that typically contains many stars, i.e., in each of the N_cells cells in which we have divided the observed Galactic region. In this way, in each cell we have the estimation of the average values of the distance, velocity components, and their corresponding errors.

4.1. Velocity Asymmetries with DR3

LIM gives estimations of the three velocity components, their errors, and rms values in the N_cells cells in which we have divided the galactic region. Figure 1 shows the maps representing the projection onto the Cartesian (X, Y)-plane of the three velocity components V_R , V_ϕ , and V_Z (from top to bottom, left panels), their errors (middle panels), and their rms values (right panels). Note that the deconvolution has included all stars with Δπ < π and ∣b∣ < 10°.

Zoom In Zoom Out Reset image size

The Galactic region within such constraints was divided into 36 line-of-sight cells, each of them with Δℓ = 10°. Then, in each of these cells we have applied the deconvolution technique discussed above. Hereafter only the cells where the number of stars is N_stars ≥ 6 are plotted. The error ${\rm{\Delta }}\overline{r}$ in the estimation of the distance $\overline{r}$ is

$$\begin{eqnarray}&&{\rm{\Delta }}\bar{r}={\sigma }_{\bar{r}},\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{\overline{r}}$ is the dispersion of $\overline{r}$ obtained by LIM. Note that this is a systematic error and it cannot be reduced by increasing N_stars in each cell.

Errors on the velocity components simply are

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{V}_{r} & = & \displaystyle \frac{\sigma ({V}_{r})}{\sqrt{{N}_{\mathrm{stars}}}}\\ {\rm{\Delta }}{\mu }_{{\ell }} & = & \displaystyle \frac{\sigma ({\mu }_{{\ell }})}{\sqrt{{N}_{\mathrm{stars}}}}\\ {\rm{\Delta }}{\mu }_{b} & = & \displaystyle \frac{\sigma ({\mu }_{b})}{\sqrt{{N}_{\mathrm{stars}}}},\end{array}\end{eqnarray} \tag{ 2 }$$

i.e., they can be reduced by increasing N_stars.

As discussed in López-Corredoira & Sylos Labini (2019), we neglect the covariance terms in the errors on μ_ℓ , μ_b , and $\overline{r}$, that is, we assume that these errors are independent of each other. Note that rms values plotted in the right panels of Figure 1 were corrected (subtracted quadratically) from the measurement errors of V_r , μ_ℓ , and μ_b . Uncertainties on V_R and V_ϕ are smaller toward the anticenter because the separation of both components is independent of the distance. Moreover, V_R only depends on V_r , so it is insensitive to the distance errors, which instead affect the determinations of both V_ϕ and V_Z .

We have tested the impact of the zero-point correction given by Lindegren et al. (2021), using the publicly available Python package, ⁶ which calculates the zero-point as a function of ecliptic latitude, magnitude, and color. By comparing Figure 2, which includes the small parallax zero-point bias, with Figure 1, we can conclude that, in line with similar tests presented in López-Corredoira & Sylos Labini (2019), such a correction has very minor effects on the derived maps in the region in the anticenter direction, where the detected measurement errors are the smallest ones, while significant differences are found where the errors are larger. However, the region toward the anticenter is the one relevant for the analysis of the velocity profiles that we present below. Indeed, one may note in the middle panels of Figures 1–2 that the error distribution shows a "horn"-like shape, delimiting the region where errors are the smallest ones: the velocity dispersion displays its lower value along X for small (in absolute value) Y, i.e., the direction of the anticenter.

Zoom In Zoom Out Reset image size

The azimuthal velocity map (see second row, left panel of Figure 1) displays patterns such that V_ϕ decreases with distance in the range [8, 25] kpc, while it increases for both Y > 0 and Y < 0. The vertical motion map (see the bottom left panel of Figure 1) presents an "arc shape" similar to the Gaia DR2 map shown in López-Corredoira & Sylos Labini (2019): this is now much better resolved, a fact supporting the robustness of LIM. Similarly, both the behavior of the error and the behavior of the dispersion (middle and right panels of Figure 1) show the same patterns as in the maps obtained from the analysis of DR2. Figure 3 shows a zoomed-in version of the radial and vertical velocity maps presented in Figure 1.

Zoom In Zoom Out Reset image size

As seen in Figure 4, it shows the edge-on projection, i.e., onto the (X, Z)-plane, of the three velocity components (as in Figure 1), with the constraint 160° < ℓ < 200°; this corresponds to the galactic anticenter region. One may note that there are several asymmetries with positive and negative velocity gradients of amplitude from 10 to 25 kpc. The errors become larger between Z = 2 and 4 kpc and between −4 and −2 kpc, and X ≤ 15 kpc. Note that in Figure 4 uncertainties are larger toward the Galactic pole because the binning with constant Δb reduces the number of sources in that region.

Zoom In Zoom Out Reset image size

Figure 5 shows the projection on the (R, Z)-plane, where the vertical coordinate has a broader range, i.e., Z ∈ [−15, 15]; one may note asymmetrical kinematics patterns such as radial non-null motions, the "horn-like" shape in azimuthal velocity, and nonzero vertical bulk motions.

Zoom In Zoom Out Reset image size

Figure 6 shows the radial velocity profile along the azimuth, in different radial bins and in the range 12−24 kpc. At least in the range between −20° and 20°, the larger is the distance, the smaller is the radial velocity; beyond this range in azimuth errors are too large to make a reliable estimation of the velocity. Figure 7 shows the vertical velocity profile along the azimuth in different radial bins and in the range 4−20 kpc. In addition, V_Z (ϕ) increases from 4 to 16 kpc, i.e., the larger is the distance, the larger is the average vertical velocity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the velocity maps and profiles (see below) obtained with the Gaia DR3 are consistent with those obtained by López-Corredoira & Sylos Labini (2019) with Gaia DR2, a fact that shows that LIM is a robust and reliable technique. This same result can be inferred by a simple visual comparison of Figure 1 with Figure 8 of López-Corredoira & Sylos Labini (2019) and of Figure 4 with their Figure 9: the maps derived from DR3 well agree with those of DR2 even when we consider the region of DR2 where errors are larger.

4.2. Asymmetric Motions: a Qualitative View

4.3. Rotation Curve and Radial Velocity Profiles

We now discuss the determination of the radial profiles of the three velocity components at different vertical heights. We then consider the derivation of the velocity rotation curve from the Jeans equation, stressing the underlying hypotheses, and then we describe in detail our results, comparing them, in particular, with those of Eilers et al. (2019).

4.3.1. Estimation of Velocity Moments

As mentioned above, LIM gives estimations of the velocity components and their dispersion, i.e., $({V}_{\alpha }^{i},\,{\sigma }_{{V}_{\alpha }}^{i})$ with α = R, ϕ, Z in the i = 1,.., N_cells cells into which we have divided the galactic region. Thus, we can estimate the average velocity components as

$$\begin{eqnarray}&&\overline{{V}_{\alpha }}=\displaystyle \frac{{\sum }_{i}\tfrac{{V}_{\alpha }^{i}}{{({\sigma }_{{V}_{\alpha }}^{i})}^{2}}}{{\sum }_{i}\tfrac{1}{{({\sigma }_{{V}_{\alpha }}^{i})}^{2}}}.\end{eqnarray} \tag{ 3 }$$

Note that we compute in bins of size ΔR = 1 kpc and at different heights in ${Z}_{\min },{Z}_{\max }$: for each R the sums in Equation (3) include all the cells that satisfy such constraints. The estimation of the variance of $\overline{{V}_{X}}$ is

$$\begin{eqnarray}&&{\sigma }_{\overline{{V}_{\alpha }}}^{2}=\displaystyle \frac{1}{{\sum }_{i}\tfrac{1}{{({\sigma }_{{V}_{\alpha }}^{i})}^{2}}},\end{eqnarray} \tag{ 4 }$$

with the same constraints as before.

4.3.2. Results for Velocity Profile

Figure 8 shows the transversal, radial, and vertical velocity profiles along the radial distance and computed in different vertical slices of size ΔZ (left panels) and in bins of size ΔZ = 1 centered at Z_c plus bins of size ΔZ = 1 centered at −Z_c (right panels). Hereafter we will refer to these two determinations as integral and differential velocity profiles.

Zoom In Zoom Out Reset image size

The average azimuthal velocity profile (Figure 8, top left panel) shows a clear monotonic decreasing trend from ∼220 km⁻¹ at R ≈ 10 kpc to ∼160 km⁻¹ at R ≈ 30 kpc. The differential determination of the azimuthal velocity displays a significant trend at small radii, i.e., for R < 15 kpc V_ϕ decreases as ∣Z_c ∣ increases (top right panel of Figure 8). Instead, for R > 15 kpc, V_ϕ does not show significant variations with Z. It should be noticed that for ∣Z∣ > 2 errors in the determination of the azimuthal velocity are large (up to ≈100 km⁻¹—see rms values in the right panel of Figure 5), and thus it might be required to confirm these results by a further data release where the error on the parallax will be lowered.

The radial velocity profile on the Galactic plane (i.e., Z_c = 0; Figure 8, middle panels) decreases as a function of radius for R < 20 kpc. Instead, the larger is Z_c , the larger is the positive value of V_R , which reaches ≈30 km⁻¹ in the outermost regions explored.

Finally, the average vertical velocity (Figure 8, bottom panels) shows a decreasing trend on the Galactic plane (i.e., Z_c = 0) for R > 15 kpc with a gradient of 20 km⁻¹ until 30 kpc.

The effect of the zero-point correction on the velocity profiles is shown in Figure 9. As noticed above, the zero-point correction affects the behavior at large radii by introducing larger errors, and for this reason such a correction reduces the range of radii where we can reconstruct kinematic properties. We find that the velocity profiles, when computed in the same galactic region, nicely overlap.

Zoom In Zoom Out Reset image size

4.3.3. The Jeans Equation

A basic assumption often used to interpret Galactic dynamics is that the disk is in equilibrium or that the gravitational potential is stationary; this hypothesis is definitely challenged by the rich complexity of velocity substructures revealed by the 3D kinematic maps provided by the Gaia mission. Theoretically it is not evident how to take into account such streaming motions in all velocity components to construct a self-consistent description of the galaxy. In what follows we will use the time-independent Jeans equation in an axisymmetric gravitational potential to compute the rotation curve. Chrobáková et al. (2020) have shown that, as long the amplitude of the radial velocity component is small compared to that of the azimuthal one, the Jeans equation provides a reasonable approximation to the system. From our analysis we may conclude that up to 30 kpc on the Galactic plane perturbations in the radial velocity should be small, and thus we may use the Jeans equation to compare observations with theoretical models.

Assuming an axisymmetric gravitational potential of the Milky Way, we use the Jeans equation (in cylindrical coordinates R, Z, ϕ; Binney & Tremaine 2008) to link the moments of the velocity distribution and the density of a collective of stars to the gravitational potential, i.e.,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \nu \langle {V}_{R}\rangle }{\partial t} & + & \displaystyle \frac{\partial \nu \langle {V}_{R}^{2}\rangle }{\partial R}+\displaystyle \frac{\partial \nu \langle {V}_{R}{V}_{Z}\rangle }{\partial z}\\ & + & \nu \left(\displaystyle \frac{\langle {V}_{R}^{2}\rangle -\langle {V}_{\phi }^{2}\rangle }{R}+\displaystyle \frac{\partial {\rm{\Phi }}}{\partial R}\right)=0,\end{array}\end{eqnarray} \tag{ 5 }$$

where ν denotes the density distribution.

The circular velocity curve in an axisymmetric gravitational potential Φ of a disk galaxy is defined as

$$\begin{eqnarray}&&{V}_{c}^{2}(R)=R\displaystyle \frac{\partial {\rm{\Phi }}}{\partial R}{| }_{z\approx 0}.\end{eqnarray} \tag{ 6 }$$

By assuming a steady state, the time-dependent term in Equation (5) is set to zero and we get

$$\begin{eqnarray}&&{V}_{c}^{2}=\langle {V}_{\phi }^{2}\rangle -\langle {V}_{R}^{2}\rangle \left(1+\displaystyle \frac{\partial \mathrm{ln}\nu }{\partial \mathrm{ln}R}+\displaystyle \frac{\partial \mathrm{ln}\langle {V}_{R}^{2}\rangle }{\partial \mathrm{ln}R}\right)+\displaystyle \frac{\partial \nu \langle {V}_{R}{V}_{z}\rangle }{\partial z}.\end{eqnarray} \tag{ 7 }$$

We assume that the volume density can be written as

$$\begin{eqnarray}&&\nu (R,z)={\rho }_{0}\exp \left(-\displaystyle \frac{R}{{h}_{R}}\right)\exp \left(-\displaystyle \frac{| z| }{{h}_{z}}\right),\end{eqnarray} \tag{ 8 }$$

where h_R is the scale length of the disk and h_z is the scale height (we use the same values of Chrobáková et al. 2020). By defining for the three velocity components

$$\begin{eqnarray}&&\langle {V}_{X}^{2}\rangle =\langle {V}_{X}{\rangle }^{2}+{\sigma }_{\langle {V}_{X}\rangle }^{2},\end{eqnarray} \tag{ 9 }$$

we find that Equation (7) can be written as

$$\begin{eqnarray}\begin{array}{rcl}{V}_{c}^{2} & = & \langle {V}_{\phi }{\rangle }^{2}+{\sigma }_{\langle {V}_{\phi }\rangle }^{2}+\left(\langle {V}_{R}{\rangle }^{2}+{\sigma }_{\langle {V}_{R}\rangle }^{2}\right)\displaystyle \frac{R-{h}_{R}}{{h}_{R}}\\ & & -2R\langle {V}_{R}\rangle \displaystyle \frac{\partial \langle {V}_{R}\rangle }{\partial R}-R\displaystyle \frac{\partial {\sigma }_{\langle {V}_{R}\rangle }^{2}}{\partial R}\\ & & +\displaystyle \frac{R}{{h}_{z}}\displaystyle \frac{z}{| z| }\langle {V}_{R}{V}_{z}\rangle -R\displaystyle \frac{\partial \langle {V}_{R}{V}_{z}\rangle }{\partial z}.\end{array}\end{eqnarray} \tag{ 10 }$$

4.3.4. Rotation Curve

Following Eilers et al. (2019) we have neglected the terms with V_Z in Equation (10), because the cross-term 〈V_R V_Z 〉 and its vertical gradient are ≈2–3 orders of magnitude smaller compared to the remaining terms; hence, their effects are negligible.

The resulting integral and differential rotation curve are reported in Figure 10 (the values for ∣Z∣ < 2 kpc are reported in Table 1). Note that we limited the analysis to the Galactic plane, ∣Z∣ < 3 kpc, i.e., in the region where the disk is highly dominant over the stellar halo contribution.

Zoom In Zoom Out Reset image size

Table 1. Measurements of the Circular Velocity of the Milky Way for the Gaia DR3 Sample in the Anticenter Region, i.e., with the Constraints 160° < ℓ < 200° and ∣Z∣ < 3 kpc

R (kpc)	vc (km s−1)	${\sigma }_{{v}_{c}}$ (km s−1)
9.5	221.3	2.5
10.5	222.6	2.4
11.5	220.5	2.4
12.5	222.9	2.4
13.5	224.1	2.2
14.5	220.7	2.3
15.5	218.1	2.3
16.5	215.5	2.4
17.5	213.0	2.6
18.5	209.4	2.9
19.5	205.4	3.1
20.5	201.8	2.7
21.5	198.4	3.3
22.5	194.3	3.5
23.5	193.7	2.4
24.5	194.4	2.1
25.5	178.5	7.4
26.5	175.5	5.8
27.5	175.3	7.7

Note. Columns show the Galactocentric radius, the circular velocity, and its error bar.

Download table as:  ASCIITypeset image

By comparing V_C with the azimuthal velocity profile (top panel of Figure 8), we notice that the additional terms from the Jeans equation only contribute as small perturbations. In particular, the main features that we have observed for the azimuthal velocity profile are present also for the rotation curve, namely, the decrease in amplitude passing from 10 to 30 kpc and, when we consider the differential measurement, the trend of decreasing amplitude at small radii with the increase (in absolute value) of the vertical height. As mentioned above (see Figure 9), we have tested that the systematic contribution due to the zero-point correction does not change significantly the behaviors observed.

That the rotation curve of the Milky Way was decreasing for R > 12 kpc was also found by several authors with different tracers and modeling (see, e.g., Dias & Lépine 2005; Xue et al. 2008; Bovy et al. 2012; Kafle et al. 2012; Gibbons et al. 2014; Lopez-Corredoira 2014; Reid et al. 2014; Galazutdinov et al. 2015; Gaia Collaboration et al. 2018; Eilers et al. 2019; López-Corredoira & Sylos Labini 2019; Mróz et al. 2019; Jiao et al. 2021; Bird et al. 2022; Wang et al. 2022). In particular, Eilers et al. (2019) considered a sample with the 6D phase-space coordinates of 23,000 luminous red giant stars, with precise parallaxes determined by combining spectral data from APOGEE DR14 with photometric information from the Wide-field Infrared Survey Explorer, 2MASS, and Gaia DR2. They measured that the circular velocity curve shows a gentle but significant decline with increasing radius and can be well approximated by a linear function up to 25 kpc:

$$\begin{eqnarray}&&{V}_{C}(R)=V({R}_{\odot })+\beta (R-{R}_{\odot }),\end{eqnarray} \tag{ 11 }$$

where R_⊙ is the distance of the Sun from the Galactic center, V(R_⊙) = 229 ± 0.2 km s⁻¹, and the slope was found to be β = −(1.7 ± 0.1) km s⁻¹ kpc⁻¹. In Figure 8 we report the rotation curve determined by Eilers et al. (2019), which nicely agrees with our estimation; in turn, their results are in good agreement with previous determinations by Kafle et al. (2012), Lopez-Corredoira (2014), and Huang et al. (2016), although these three determinations have larger error bars (see Figure 3 of Eilers et al. 2019).

The result of Eilers et al. (2019) is in reasonably good agreement with another recent analysis of the Milky Way's circular velocity curve by Mróz et al. (2019) (also reported in Figure 8 for comparison). This was based on a sample of 773 classical Cepheids with precise distances based on mid-infrared period−luminosity relations coupled with proper motions and radial velocities from Gaia. They found, in the range of radii between 5 and 20 kpc, a somewhat larger slope of β = −(1.4 ± 0.1) km s⁻¹ kpc⁻¹. Note from Figure 8 that, however, the number of Cepheids for R > 15 kpc drops significantly, so that the slope β was measured on a very limited range of radii.

By extending the rotation curve to 27.5 kpc, we find that β = −(2.3 ± 0.2) km s⁻¹ kpc⁻¹, which is smaller than the other two determinations mentioned above. However, in the range of radii where they overlap, i.e., for R < 20 kpc for the determination by Mróz et al. (2019) and for R < 25 kpc for the one by Eilers et al. (2019), the three measurements are in reasonably good agreement with each other. In this respect we note that our result is independent of those of Mróz et al. (2019) and Eilers et al. (2019), as we used a different data set and a different method to determine the rotation curve, namely, we have constructed a coarse-grained sample and have applied to it a statistical method to measure the average velocity components and their dispersion, which was tested to provide reliable results, rather than analyzing the velocities of individual stars.

In order to estimate systematic uncertainties on the circular velocity curve arising from our data sample, we split the galactic region into two disjoint smaller portions, one with 0 kpc < Z < 2 kpc and the other one with −2 kpc < Z < 0 kpc (a similar result is obtained by dividing the sample into b > 0° and b < 0°). We then computed the rotation curve in the two disjointed regions and estimated the systematic uncertainties on the circular velocity by the difference between the resulting fit parameters from the two disjoint data sets. We find that the systematic error on the slope of the rotation curve β (see Equation (11)) is ≈0.4 km s⁻¹ kpc⁻¹, which corresponds to an uncertainty of the order of 20%. An additional contribution to the systematic comes from the fact that, for calculating the circular velocity curve from Equation (7), we neglected the cross terms 〈V_z V_R 〉. The large noise affecting this term does not allow a precise determination of its effect, but in agreement with Eilers et al. (2019), we estimate it to be of the order of a few percent.

4.3.5. Discussion

The two determinations of the Galaxy rotation curve by Mróz et al. (2019) and Eilers et al. (2019) and the one presented in this work (which up to 20 kpc coincides with the one by López-Corredoira & Sylos Labini 2019) are different from others reported in the literature (see, e.g., Bhattacharjee et al. 2014; Sofue 2020), where the rotation curve of the Milky Way did not present a decrease in the range of distances between 15 and 30 kpc. In such estimations V_C (R) was measured when, except in a few cases, the full 3D velocity information of the tracers was not available, and it has to be reconstructed from only the measured line-of-sight velocity and positional information of various tracer objects in the Galaxy. On the other hand, results by Mróz et al. (2019), by Eilers et al. (2019), and by us in this work share the key fact that they are based on different measurements that have an unprecedented precision and accuracy in the determination of the distances and that allow us to have information on the 6D phase space. The knowledge of the phase-space distribution allows us to control for possible systematic effects of a different kind, such as that induced by coherent radial and vertical motions. Given the agreement between these determinations, we conclude that our results are reliable and that possible systematic effects in the determination of the rotation curve should be marginal; of course, when the next data release of the Gaia mission is published, it will be possible to further test the rotation curve for R > 20 kpc. This situation actually shows that the better accuracy of the stars' distances, together with the great amount of new information provided by the Gaia mission, is the key improvement to our knowledge of the Galaxy, its kinematics, and its dynamics.

Concerning the estimation of the mass of the Milky Way, we note that Eilers et al. (2019), by considering the standard Navarro−Frenk−White halo model, found a virial mass of M = (7.25 ± 0.25) × 10¹¹ M_⊙, which is significantly lower than what several previous studies suggest. The declining trend of the rotation curve continues when we extend the range of radii from 25 to 27.5 kpc, which implies that the best fit with a Navarro−Frenk−White halo model (Navarro et al. 1997) should correspond to an even smaller Milky Way virial mass (less massive than Gaia EDR3 work by about 20%). A detailed discussion of the dynamical implication of these results will be presented in a forthcoming publication.

Another noticeable feature of the circular velocity that we have detected, namely, the trend of the decreasing amplitude for R < 15 kpc with the increase (in absolute value) of the vertical height, has, to our knowledge, not been noticed before. A similar study was presented in Chrobáková et al. (2020) with Gaia DR2, reaching distance up to R ≈ 20 kpc and heights only up to ∣Z∣ < 2 kpc, where these trends are less noticeable. Moreover, the binning of data in Chrobáková et al. (2020) was finer, meaning that the derivatives in Equation (10) are more influenced by fluctuations, making them less reliable. Thus, the rotation curves of Chrobáková et al. (2020) are less robust than our current analysis, and the trends we see now were not noticed then. As for the case of the determination of V_ϕ (R, Z), we emphasize that for ∣Z∣ > 2 kpc errors in the determination of the azimuthal velocity are large (up to ≈100 km⁻¹—see rms values in the right panel of Figure 5).

We have used LIM to analyze the Gaia DR3 data set (Gaia Collaboration et al. 2022b). LIM can solve the deconvolution of large Gaussian errors that affect the measurements of stellar distance, and it was previously applied to the DR2 data set (López-Corredoira & Sylos Labini 2019); in this way it was possible to derive the kinematic maps of the Galaxy covering a region where the relative error in distance was larger than 20%, thus extending the range of distances for the kinematic analyses with respect to those presented by Gaia Collaboration et al. (2018). The new analysis discussed in the present work allows us to explore a range of Galactocentric distances up to ≈30 kpc, while the range of distances covered by studies using only stars with distance errors <20% also increased, by passing from DR2 to DR3, from 13 (Gaia Collaboration et al. 2018) to 18 kpc (maximum value) (Gaia Collaboration et al. 2022a).

The first noticeable result is that LIM applied to the DR3 data is compatible with the results of DR2 (López-Corredoira & Sylos Labini 2019). This is already an important result showing that LIM converges, that is, by lowering the parallax errors and by increasing the number of sources, i.e., the measured stars. The results well confirm those obtained earlier in a data set with fewer objects and larger errors. As the method is designed to work when the distribution of error is Gaussian, this means that the observational parallax errors satisfy such a condition.

In addition, as a second key result of our work, we find that the new extended maps of the Galactic disk cover the regions of the outer disk that are farther from the Galactic center, whose stars reach R ≈ 30 kpc (López-Corredoira et al. 2018). These maps show that there are large-amplitude and coherent streaming motions in all velocity components. In particular, the radial velocity profile shows an increase toward the outermost region of the Galaxy, but off plane, with a detected value of V_R ≈ 30 km s⁻¹. The azimuthal velocity also shows a clear decreasing trend as a function of radius. In addition, we have found a marked change of V_ϕ (R) at small radii, i.e., R < 15 kpc, when we consider its determination at different heights. In summary, the new extended maps confirm that the Galaxy kinematics is characterized by significant coherent streaming motions in all velocity components as found, at smaller radii, by, e.g., Gaia Collaboration et al. (2018), Antoja et al. (2018), López-Corredoira & Sylos Labini (2019), and Khoperskov et al. (2020).

By computing the rotation curve through the Jeans equation, assuming that the Galaxy is in a steady state and that the galactic potential is axisymmetric, we found that V_c (R)/V_ϕ (R) shows a significant decline from ≈15 to 30 kpc of more than 50 km s ⁻¹. This result is in reasonable agreement with the recent findings by Mróz et al. (2019) up to 20 kpc and by Eilers et al. (2019) up to 25 kpc and extends them to 27.5 kpc but with smaller errors. (Of course, the behavior of the rotation curve in this paper agrees with that found by López-Corredoira & Sylos Labini 2019 by applying LIM to the Gaia DR2 sources.) These three results used different samples that have an unprecedented precision and accuracy in the determination of the distances and that allow us to have information on the 6D phase space. Another interesting result that we found is that the rotation curve, as well the azimuthal velocity, presents a marked dependence on the height for R < 15 kpc, whereas, in both cases, at larger R the dependence on Z is negligible. In a forthcoming work we will interpret these behaviors in light of different galactic models providing an estimation of the Galaxy's mass.

We would like to thank Roberto Capuzzo-Dolcetta very much for useful comments and insightful discussion. H.-F.W. acknowledges the support from the project "Complexity in self-gravitating systems" of the Enrico Fermi Research Center (Rome, Italy) and science research grants from the China Manned Space Project, Nos. CMS-CSST-2021-B03 and CMS-CSST-2021-A08. Z.C. was supported by VEGA—the Slovak Grant Agency for Science, grant No. 1/0761/21, and by the Erasmus+ programme of the European Union under grant No. 2020-1-CZ01-KA203-078200. M.L.C. were supported by the grant PGC-2018-102249-B-100 of the Spanish Ministry of Economy and Competitiveness (MINECO). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.