To provide a cleaner sample for further analysis we experimented with the
data, reaching the conclusion that retaining only DCEPs with distances smaller
than 30 kpc and relative distance uncertainty better than 10% is a good
compromise between precision and completeness. This selection removed 240
DCEPs from the Gaia DR3 sample and 230 literature DCEPs, leaving us with a
total final dataset of 3306 DCEPs useful for further analysis. The
$G$-magnitude distribution of this selected sample is shown in Fig. 2. The
number of faint objects, especially among the literature sample is now
drastically reduced.

Figure 2: Histogram of the $G$ magnitude for the DCEP sample used in this
work. We plot in cyan and red the histograms for the objects from the Gaia DR3
catalogue and those taken from the literature, respectively. The thick line
shows the magnitude distribution of the DCEPs after the selection in age
($<$200 Gyr), while the thin line shows the histogram of the objects selected
in age that also possess a line-of-sight measurement. In the last two cases
the Gaia DR3 and literature samples were merged.

####  2.2.2 Line-of-sight velocities.

The spectra collected with the RVS spectrometer on board Gaia allows the
measurement of time series line-of-sight velocity (RV) values for millions of
stars with $G$ magnitude brighter than $\sim$15-16 mag (for details see sect.
6.4.8 of Sartoretti et al., 2022). In DR3 the RV time series are released for
a selected sample of 774 DCEPs666Among these objects, 15 and 9 belong to the
Large and Small Megallanic Clouds, respectively. For this subsample, the SOS
Cep&RRL pipeline computed average RV values and relative uncertainties by
fitting the RV curves folded according to the stars’ periods and makes them
available in the vari_cepheid table (see Clementini et al., 2022; Ripepi et
al., 2022b, for details). Nevertheless, for a much larger number of DCEPs in
our sample, mean RVs estimated from the spectroscopic pipeline are available
in the main gaia_source table. We checked whether these arithmetically
computed mean RV values are usable for variable objects such as DCEPs by
comparing the RV values in the vari_cepheid and the gaia_source tables. We
found a perfect agreement with a mean difference of 0.6$\pm$6 km/s and no
visible trend (see also Clementini et al., 2022). On this basis we decided to
use the gaia_source catalogue RVs for our DCEP sample. In total we have RV
estimates for 2059 DR3 and 67 literature DCEPs, respectively. The
uncertainties on these values can be evaluated on the basis of fig. 6.13 by
Sartoretti et al. (2022).

####  2.2.3 Ages of DCEPs

It has been known for a long time that the DCEPs follow a period-age (PA)
relation (see e.g. Bono et al., 2005; Anderson et al., 2016, and references
therein). More recently, De Somma et al. (2021), based on an updated
theoretical pulsation scenario, devised a more accurate period-age-metallicity
(PAZ) relation. Since the DCEP PWZ relation allows us to obtain individual
accurate distances, it follows that the DCEP PAZ allows us to date whatever
region of the Galactic disk where a DCEP is present. To take advantage of this
powerful tool, we adopted the following equation (see Table 9 of De Somma et
al., 2021):

$\displaystyle\log t=(8.423\pm 0.006)-(0.642\pm 0.004)\log P$ (3)
$\displaystyle-(0.067\pm 0.006){\rm[Fe/H]}$

with rms=0.081; $t$, $P$ and [Fe/H] are the age (years), the period (days),
and the iron abundance (dex). This relation is valid for F-mode pulsators and
it was calculated using evolutionary tracks including overshooting. As the PAZ
relation is not available for 1O-mode DCEPs (see De Somma et al., 2021, for
full details), we decided to fundamentalise their periods according to the
Feast & Catchpole (1997) relation: $P_{F}=P_{1O}/(0.716-0.027\log P_{1O})$,
where $P_{F}$ and $P_{1O}$ are the F and 1O mode DCEP periods, respectively.
In this way we were able to calculate the ages for every DCEP in our sample.
As we want to use the DCEPs to trace the MW arms, we decided to use only DCEPs
younger than 200 Myr. Therefore the sample used in the following is composed
of 2808 pulsators, 1948 of which also have line-of-sight velocity
measurements. Table 2 shows selected DCEPs properties which are not published
in other Gaia catalogues or papers.

Table 2:  Selected parameters for the DCEPs sample.

source_id | $\ell$ | $b$ | $d$ | $\mu$ | $\sigma_{\mu}$ | [Fe/H] | $\sigma_{\mathrm{[Fe/H]}}$ | flag | source | logAge | $\sigma_{\mathrm{logAge}}$
---|---|---|---|---|---|---|---|---|---|---|---
| (deg) | (deg) | (kpc) | (mag) | (mag) | (dex) | (dex) |  |  | (dex) | (dex)
4060910068247394432 | 0.38444 | 1.96190 | 8.02 | 14.522 | 0.129 | 0.50 | 0.20 | 1 | P21 | 8.594 | 0.016
4049125051634137600 | 0.45968 |  $-$5.42326 | 14.46 | 15.801 | 0.117 | 0.19 | 0.20 | 1 | P21 | 8.706 | 0.015
4048895253682114432 | 0.47058 |  $-$6.65562 | 19.14 | 16.410 | 0.114 |  $-$0.05 | 0.20 | 1 | Gaia_DR3 | 8.623 | 0.015
4056461478623363968 | 359.98638 |  $-$1.45395 | 15.81 | 15.995 | 0.157 | 0.11 | 0.20 | 1 | Gaia_DR3 | 7.715 | 0.015

777The full table is available as an electronic table via CDS. The source_id
is the Gaia source_id, $\ell$ and $b$ are the galactic coordinates; $d$ is the
heliocentric distance in kpc; $\mu$ and $\sigma_{\mu}$ are the distance
modulus and its uncertainty; [Fe/H] and $\sigma$[Fe/H] are the iron abundance
and its uncertainty; flag=0 or 1 means that the metallicity was taken from the
astrophysical parameters or calculated from the metallicity gradient of the
Galactic disc; source lists the provenance of the DCEP source: Gaia_DR3 means
that the star is included in the Gaia DR3 vari_cepheids catalogue, P21 or Inno
indicate that the DCEPs were taken from the Pietrukowicz et al. (2021) or Inno
et al. (2021) catalogues, respectively; logAge and $\sigma_{\mathrm{logAge}}$
are the decimal logarithm of the age and its uncertainty.

###  2.3 OB stars

To select young stars on the upper main sequence we use the effective
temperatures provided in DR3, selecting for stars with
${\mathrm{T_{eff}}}>10\,000$K. For hot stars there are in general two sets of
effective temperatures provided in DR3. One set is provided by a General
Stellar Parameterizer from Photometry (GSP-Phot, see Andrae et al., 2022),
that estimates stellar parameters using the Gaia
$G_{\mathrm{BP}}$/$G_{\mathrm{RP}}$ spectrophotometry, astrometry and $G$ band
photometry. GSP-Phot makes different sets of parameter estimates using
different stellar libraries, then for each source chooses one of these as the
”best”. Here we use this set of ”best” parameters, as reported in the main
Gaia source table. Another set of parameters is estimated from a software
module (ESP-HS) that was optimised specifically for hot stars and which uses
the BP/RP spectrophotometry, without the astrometry, together with the RVS
spectra if also available (Creevey et al., 2022; Fouesneau et al., 2022). This
second set of parameters is made available in the astrophysical_parameters
table. Because of different quality filters for these different methods, Gaia
sources may have one or both sets of effective temperatures, or remain without
a temperature estimate. Indeed, only about half the sample of stars with
${\mathrm{T_{eff}}}>10\,000$K from either method have temperatures from both.
From a detailed comparison of those sources with stellar parameters from one
or both methods, we settled on the following criteria:

  * •

For the stars with only GSP-Phot temperatures we use the spectral type
determined by ESP-HS for all sources with Gaia BP/RP spectrophotometry as an
additional assurance of quality. That is ${\mathrm{T_{eff}}}>10\,000$K and the
ESP-HS spectral type flag set to ”O”, ”B” or ”A”.

  * •

For the stars with only ESP-HS temperatures we require that the effective
temperature be in the range $10\,000<{\mathrm{T_{eff}}}<50\,000$K, as it was
found that the small fraction of sources with ${\mathrm{T_{eff}}}>50\,000$K
are likely to be unreliable (Fouesneau et al., 2022).

  * •

For sources with both sets of stellar parameters we require that the effective
temperature ${\mathrm{T_{eff}}}>8000$K for GSP-Phot and
${\mathrm{T_{eff}}}>10\,000$K for ESP-HS, letting the confirmation from GSP-
Phot verify the sources with ESP-HS hotter than 50 000K.

We note that we only use the effective temperature for selection. A comparison
of our temperatures against those found in the literature (Mathur et al.,
2017; Abolfathi et al., 2018; Xiang et al., 2021) for stars in the range 8000
to 10 000K show an RMS difference of less than 900K and offsets less than
400K. While the differences increase for higher ${\mathrm{T_{eff}}}$ measured,
they remain small enough to that we remain confident that they should
nevertheless be in our sample. We also note that the use of ESP-HS products
for all three of the cases above effectively poses an apparent magnitude limit
on this sample of 17.65 in $G$.

As a temperature selection will introduce undesired sub-dwarfs and white
dwarfs in our sample, we impose the additional criterion
$G+5\log(\varpi/100)<2.+1.8({G_{\mathrm{BP}}}-{G_{\mathrm{RP}}})$ to remove
sources fainter than our target upper main sequence stars. The colour term
takes into account extinction, where 1.8 is approximately the slope of the
reddening vector in $(G,G_{\mathrm{BP}}-G_{\mathrm{RP}})$ space. To capture
distant sources with negative parallaxes we rewrite this criteria in the form

$(\varpi/100.)^{5}<10.^{(2.-G+1.8*({G_{\mathrm{BP}}}-{G_{\mathrm{RP}}}))}\,.$
(4)

(See appendix A for an example query.) Together the above criteria give us 923
700 stars, but we find that outside the plane of the Galaxy we have a
significant number of stars in the direction of the Large (LMC) and Small
(SMC) Megallanic Clouds, as well as a number of globular cluster members. We
remove these contaminants by keeping only stars whose distance from the
galactic plane is less than 300pc. This reduces the sample to 621 609 stars.
Finally we use the astrometric fidelity indicator $f_{a}$ (with values
$0<f_{a}<1$) of Rybizki et al. (2022), to remove sources with suspect
astrometry, keeping stars with $f_{a}>0.5$ as recommended in Zari et al.
(2021). This last criteria removes only 7% of the sample, leaving us with a
final sample of 579 577 stars, of which 91 836 (15.8%) have line-of-sight
velocities. However, we note that a fraction of these have line-of-sight
velocities estimated using an RVS spectral template with temperatures very
different from the effective temperatures finally estimated for them, and
therefore likely to have errant line-of-sight velocities (Blomme et al.,
2022). Removing stars whose rv_template_temp) is less than 7000K gives us 77
659 stars with valid line-of-sight velocities. Figure 3 shows the $G$
magnitude distribution of our OB sample.

Figure 3: $G$ magnitude distribution of the selected giants (red) and OB
(blue) samples. Black and grey histograms show the subsamples with line-of-
sight velocities, for the giants and OB stars respectively. The magenta
histogram are OB stars with line-of-sight velocities estimated with templates
having a ${\mathrm{T_{eff}}}$ ¿ 7000K, and is the selected sample for mapping
the velocities of the OB stars.

For the purpose of mapping we need distance estimates for our sources. While
Gaia provides parallaxes, about 40% of our sample have significant
($\sigma_{\varpi}/\varpi>0.20)$ parallax uncertainties, so that a simple
inversion of the parallax cannot be considered reliable (Bailer-Jones, 2015;
Luri et al., 2018). Figure 4 shows the distribution of
$\varpi/\sigma_{\varpi}$ for our sample as well as the subsample with line-of-
sight velocities. For the subset of our stars with GSP-Phot parameters, Gaia
DR3 also gives us distance estimates based on both astrometry and photometry.
However, 43% of our sample are sources that have only temperature estimates
from ESP-HS, and therefore no distance estimate from GSP-Phot. We therefore
adopt the ”photogeo” distances from Bailer-Jones et al. (2021) that are also
based on both astrometric and photometric data, as recommended in Fouesneau et
al. (2022), and which are available for our entire sample888These distances
can be found in the Gaia Archive as the external table,
external.gaiaedr3_distance..

Figure 4: Distribution of $\varpi/\sigma_{\varpi}$ for the OB sample (blue)
and the subsample with line-of-sight velocities (orange). The vertical red
line is at $\varpi/\sigma_{\varpi}=5$.

We mention that an alternative selection of high-fidelity OB stars is
presented in Gaia Collaboration et al. (2022a). The young B-stars from that
sample are used for a basic modelling of the Milky Way rotation curve which
results in parameters consistent with the mean OB star $V_{\phi}$ curve
derived below in section 5.

###  2.4 Giants

To select stars on the Red Giant Branch (RGB), we use the effective
temperatures and surface gravities provided in DR3, selecting for stars with
$3000<{\mathrm{T_{eff}}}<5\,500$K and ${\mathrm{logg}}<3.0$, as provided by
GSP-Phot (Andrae et al., 2022) in the main Gaia source table. (See appendix A
for an example of the query used.) These are given as the ”best” set of
parameters using a multi-spectral library approach, which for the RGB
correspond to either the MARCS and PHOENIX libraries (Fouesneau et al., 2022).
The Kiel diagram for these sources is shown in Fig. 5. We refer to this set as
the full RGB sample and it consists of $11\,576\,957$ sources. The magnitude
$G$ distribution for the full RGB sample is shown in Fig. 3, together with the
RGB sample with RVS line-of-sight velocities.

Figure 5: Kiel diagram for the selected RGB sample with astrometry.

As in the OB sample, and in order to perform density and velocity maps, we
need to choose a distance estimator. The distribution of
$\varpi/\sigma_{\varpi}$ for the RGB sample is very similar to that of the
full OB sample shown in Fig. 4, with about 39.5% of the giant sample having a
$\varpi/\sigma_{\varpi}<5$. All stars in the giant sample have GSP-Phot
parameters, thus, one option can be to use the provided distance in Gaia DR3
Andrae et al. (2022). Another option is the ”geo” and ”photogeo” distances
from Bailer-Jones et al. (2021, hereafter CBJ2021). In order to choose the
appropriate distance estimator for the large extent of the RGB sample, we
crossmatch the recent catalogue of red clump stars by APOGEE DR17 (Abdurro’uf
et al., 2022) with our RGB sample, using sky coordinates and a radius of $1$
arcsec, resulting in a common sample of $18\,322$. By comparing the absolute
difference between the reported photometric distance and the three different
possibilities mentioned above of the stars in common, we observe that the
distance estimator with less bias and dispersion is the ”photogeo” distance by
Bailer-Jones et al. (2021) (see Fig. 6). Also, and as shown in Babusiaux et
al. (2022), the large parallax uncertainties and the prior used makes the
derived GSP-Phot distances to concentrate in density forming a ring around $2$
kpc, making them inappropriate to study the inner disc. Therefore, and as in
the OB sample, we adopted the ”photogeo” as distance estimator for the RGB
sample.

Figure 6: Absolute difference between the APOGEE Red Clump distance and three
distance estimators considered in this work, GSP-Phot (left), ”geo” (middle)
and ”photogeo” (right).

In addition, as in the selection of OB stars, we keep only sources with good
astrometry, that is, with an astrometric fidelity $f_{a}>0.5$. This criteria
removes only $14\%$ of the sample, leaving us with a sample of $9\,959\,807$
stars, of which $6\,586\,329$ have line-of-sight velocities.

We found that the above selection also includes many sources from the Large
and Small Magellanic Clouds, as well as a number of globular clusters. Since
the goal of the paper is to study the Galactic disc, in the following we
perform an extra cut on the altitude with respect to the Galactic plane, i.e.
$|Z|<1$ kpc, see Sect. 3.1. We have also explored whether this selection could
include sources from the Sagittarius Dwarf galaxy, which could bias our
kinematic study. We find no evidence of Sagittarius Dwarf sources in the
proper motion map, and while red clump stars in the dwarf galaxy have a
magnitude range of $17-18$ mag (e.g., Antoja et al., 2020), which falls in the
faint end of our sample, so will be a negligible fraction in our full sample,
and are not expected in the sample with RVS line-of-sight velocities, so we do
not make any attempt for removing them.

The final RGB sample consists then of $8\,727\,344$ sources, with
$5\,730\,578$ with RVS line-of-sight velocities.

In Table 3 we summarise the number of stars in each of the selected tracers
used throughout the paper. We provide both the number of sources with full 5D
astrometry and the subsample including the 6D (astrometry and line-of-sight
velocities).

Table 3: Number of stars according to the selected tracer. Sample | Astrometry | 6D
---|---|---
Clusters (all ages) | 2681 | 2162
Young clusters (age ¡ 100 Myr) | 988 | 698
Classical Cepheids (DCEP) | 3312 | 2127
Classical Cepheids (age ¡ 200 Myr) | 2812 | 1949
Young field stars (OB) | 579577 | 77659
Giant field stars (RGB) | 8727344 | 5730578

##  3 From observations to 6D phase space

###  3.1 Mapping to configuration space

To map our tracers in a 3D Cartesian coordinate space we must transform
astrometric angular measurements $(\varpi,\alpha,\delta)$ to associated
lengths. While each of these measurements have associated uncertainties, we
are here concerned with objects reaching large distances (i.e. small
parallaxes). In this case our positional uncertainties are completely
dominated by the uncertainties in parallax (whose associated positional
uncertainty is proportional to the square of the distance), and we can safely
ignore the uncertainties (and correlations) with the angular positional
measurements. Our problem is now just reduced to determining the heliocentric
distance to each tracer and its uncertainty. In the previous section the
distance estimate adopted for each tracer population is described: for the OB
and RGB samples we adopt the ”photogeo” distances of CBJ2021, for the clusters
the inverse of their median parallax, and for the Cepheids a photometric based
distance.

It has been known for a long time that some MW DCEPs are members of OCs (see
e.g. Anderson et al., 2013, and references therein). As in this paper we use
both DCEPs and OCs to study the spiral arm properties of the MW, it is useful
to compare the distances inferred for DCEPs as explained in Sect. 2.2 with
those of their host OCs, which are based on the median parallax of the OC
members. To select possible DCEPs belonging to OCs we cross-matched the DCEPs
list with that of all the known OC members. The cross-match was carried out
using Gaia identifiers and returned 25 matches. The distance comparison is
shown in Fig. 7. The overall agreement is very good, approximately below
1$\sigma$ in all the cases, except for the farthest OC of the sample, namely
UBC 608 for which the discrepancy is however less than 1.5$\sigma$.

Figure 7: Comparison of the distances calculated for DCEPs and their host OCs.

Now that we have heliocentric distances $d$ we can easily derive heliocentric
Cartesian coordinates $(x,y,z)$ under the assumption of being in a (non-
relativistic) 3D Euclidean geometry. Since we are assuming the positions as
known, we can take advantage of the provision of $(l,b)$ of each tracer in the
Gaia Archive and use the usual transformations

$\begin{pmatrix}x\\\ y\\\ z\end{pmatrix}=d\begin{pmatrix}\cos b\cos l\\\ \cos
b\sin l\\\ \sin b\end{pmatrix},$ (5)

where positive $x$ is toward the Galactic centre. galactocentric Cartesian
coordinates are then typically derived as a simple translation of the origin:

$\begin{pmatrix}X\\\ Y\\\ Z\end{pmatrix}=\begin{pmatrix}x-R_{\odot}\\\ y\\\
z+Z_{\odot}\end{pmatrix}$ (6)

where $R_{\odot}$ and $Z_{\odot}$ are the distance of the Sun from the
Galactic centre and above the Galactic ($Z=0$) midplane.

For our choice of $R_{\odot}$ we use the geometrical determination based on
the line-of-sight velocity and relative astrometry of the resolved SagA*-S2
binary, as measured by the latest contribution of the GRAVITY Collaboration
(Gravity Collaboration et al., 2022), namely $R_{\odot}=8277\pm 9$(stat) $\pm
30$(sys) pc, though we note that this value does not agree within the
uncertainties with the independent determination of $R_{\odot}$ by Do et al.
(2019) (7.959 +/- 59 (stat) +/- 32 (sys)), nor with their previous
determinations, indicating that our assumed $R_{\odot}$ may in fact be in
error by as much as 200 to 300 parsecs. Our adopted value of $R_{\odot}$
assumes of course that the position of SagA* marks the Galactic centre, which
is expected from dynamical considerations: A supermassive black hole not
already at the centre of a large stellar system will eventually migrate to the
centre due to dynamical friction (Gualandris & Merritt, 2008). Most recently
Leung et al. (2022) have independently determined $R_{\odot}=8.23\pm 0.12$kpc
based on observed stellar kinematics toward the Galactic centre, a value that
is consistent with our assumed value, but with a more realistic estimate of
its uncertainty. For mapping in configuration space, any systematic error in
$R_{\odot}$ only results in a trivial offset in the maps. However, as
discussed further in section 3.4, it can introduce rather undesirable effects
when mapping the velocities in galactocentric cylindrical coordinates.

It should be noted that the transformation to galactocentric coordinates (Eq.
6) is an approximation, as it assumes that the $b=0$ and $Z=0$ plane are
parallel to each other. While this was the original intent when the galactic
coordinate system was defined, (Gum et al., 1960; Blaauw et al., 1960), there
may well be a residual offset due to the Sun’s height above the $Z=0$ plane.
Indeed, it was already noted at the time that determinations of $Z_{\odot}$
from Hydrogen radio emission do not coincide with those based on nearby
stellar samples. In their pre-Gaia review of our knowledge of the Milky Way,
Bland-Hawthorn & Gerhard (2016) found that estimates for the distance of the
Sun above the midplane fall between 20 and 30 pc, while more recent estimates
have generally been smaller (Yao et al., 2017; Widmark, 2019; Anderson et al.,
2019; Reid et al., 2019). However, the evidence of vertical oscillations in
the disc of the Milky Way (Bennett & Bovy, 2019) and evidence of its
disequilibrium state (Antoja et al., 2018) rather complicates this discussion,
as the local stellar mid-plane of the Galaxy could very well not coincide with
a $Z=0$ plane as defined by the average vertical density distribution of the
inner disc. Indeed, different vertical modes may be present in the gas (and
star formation tracers) and in the stars, or even between different stellar
populations, explaining some of the observed variance between the different
determinations of $Z_{\odot}$. As a case in point, using a very local sample
of stars, Gaia Collaboration et al. (2021d) find that $Z_{\odot}$ varies from
-4pc to 15pc for young to older stellar populations. The observed small
negative offset of SagA* from $b=0$ also suggests that there is a residual
tilt of $\delta\theta$ between the galactic ($b=0$) plane and the Galactic
($Z=0$) plane of the order of only $0.1^{\circ}$ (see also the discussion in
Bland-Hawthorn & Gerhard, 2016). However, since for mapping the large scale
asymmetries in the disc we are primarily concerned with the positions and
velocities of our tracers in the $(X,Y)$ plane, where the effect of this tilt
is negligible, we conveniently assume $Z_{\odot}=0$, so that $Z=z$.

###  3.2 mapping to velocity space

To map the velocities we must now use the spectroscopically measured line-of-
sight velocities $v_{r}$ together with the measured proper motions and the
distance estimator described in the previous section. As mentioned above, the
OB and RGB samples with measured line-of-sight velocities are 77 659 and 5 730
578 in number, respectively. In addition, the brighter magnitude limit of the
OB sample with line-of-sight velocities also means that the area that can be
mapped by the OB stars is much smaller than that of the RGB sample. (See
figure 8).

Figure 8: Distribution of the heliocentric distances from Bailer-Jones et al.
(2021), of the RGB (red) and OB stars (blue) with line-of-sight velocities.

The relative velocity components in the heliocentric Cartesian coordinates
defined by Eq. 5 are

${\bf{v}}_{rel}=\begin{pmatrix}u\\\ v\\\
w\end{pmatrix}=A^{\prime}_{G}\,A\,\begin{pmatrix}4.74047\leavevmode\nobreak\
{\mu_{\alpha^{*}}}\leavevmode\nobreak\ d\\\ 4.74047\leavevmode\nobreak\
{\mu_{\delta}}\leavevmode\nobreak\ d\\\ v_{r}\end{pmatrix},$ (7)

where $(\mu_{\alpha^{*}},\mu_{\delta}$) are the proper motion components,
$v_{r}$ the line-of-sight (i.e. radial) velocity, $A^{\prime}_{G}$ is the
transformation matrix from equatorial to galactic coordinates, as given by
equation 4.62 of the Gaia EDR3 online documentation, and the matrix $A$ is the
normal triad at the star999Alternatively one could first convert the
astrometry to galactic coordinates, removing $A^{\prime}$ from Eq. 7 and
substituting the proper motion components in ICRS with those in Galactic
coordinates, and then computing $\bf{v}_{rel}$ from the Galactic coordinates.
In this case the triad (Eq. 8) would be in $(l,b)$ rather than
$(\alpha,\delta)$.:

$A=\begin{pmatrix}-\sin\alpha&-\sin\delta\cos\alpha&\cos\delta\cos\alpha\\\
\cos\alpha&-\sin\delta\sin\alpha&\cos\delta\sin\alpha\\\
0&\cos\delta&\sin\delta\\\ \end{pmatrix}$ (8)

We note that for the line-of-sight velocities $v_{r}$ of the OB stars, we
applied the following correction (as prescribed in Blomme et al., 2022):

$v_{\mathrm{l}os}=\href
https://geapre.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html#gaia_source-
radial_velocity-7.98+1.135\href
https://geapre.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html#gaia_source-
grvs_mag.$ (9)

This is applied to stars where $8500\leq\href
https://geapre.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html#gaia_source-
rv_template_teff\leq 14\,500$ K and $6\leq\href
https://geapre.esac.esa.int/archive/documentation/GDR3/Gaia_archive/chap_datamodel/sec_dm_main_tables/ssec_dm_gaia_source.html#gaia_source-
grvs_mag\leq 12$.

The velocities ${\bf{v}}_{rel}=(u,v,w)$ derived above are relative to the Sun.
To put them in a galactocentric reference frame we must add the Sun’s velocity
with respect to the galactic centre, $\bf{v}_{\odot}$:

${\bf{v_{*}}}={\bf{v}}_{rel}+{\bf{v}}_{\odot}\,.$ (10)

Traditionally the Sun’s velocity $\bf{v}_{\odot}$ has been estimated from the
solar motion with respect to a local standard of rest (LSR) and an adopted
value of the velocity of the LSR, usually assumed to be in circular motion
about the galactic centre. However, thanks to the recent precise measurement
of the proper motion of the SagA*, together with $R_{\odot}$, the azimuthal
components of the Sun’s galactocentric velocity can be derived in a more
direct and precise way. From Reid & Brunthaler (2020) we have
$(\mu_{l},\mu_{b})=(-6.411\pm 0.008,-0.219\pm 0.007)$ $\rm\,mas\,yr^{-1}$ for
the proper motion of SagA*, which together with $R_{\odot}$ gives the Sun’s Y
and Z-velocity components. Meanwhile, the same reduction of the SagA*-S2 data
by the GRAVITY Collaboration that yielded $R_{\odot}$ also yields the line-of-
sight velocity toward SagA*, interpreted as the reflex motion of the Sun’s
velocity toward the Galactic centre. This gives us:

${\bf{v}}_{\odot}=\begin{pmatrix}9.3\pm 1.3\\\ 251.5\pm 1.0\\\ 8.59\pm
0.28\end{pmatrix}\rm\,km\,s^{-1}$ (11)

if we assume that SagA* is stationary with respect to the galactic centre.
(see Drimmel & Poggio, 2018, for further discussion on this approach to
deriving $\bf{v}_{\odot}$). The uncertainties in ${\bf{v}}_{\odot}$, as well
as any error in our adopted $R_{\odot}$, gives us a systematic error common to
all our galactocentric velocities ${\bf{v_{*}}}$. See section 3.4 below for
further discussion.

One should note that the $(u,v,w)$ components of our galactocentric velocities
${\bf{v_{*}}}$ are rigorously in the same coordinate system defined by Eq. 5,
so slightly tilted with respect to a $(U,V,W)$ coordinate system whose $U$-$V$
plane is parallel to the mean $Z=0$ plane of the Galaxy by the angle
$\delta\theta$ mentioned above, if $Z_{\odot}\neq 0$. However, the systematic
error introduced by ignoring this (unknown) tilt is much smaller than the
systematic error introduced by the uncertainties in ${\bf{v}}_{\odot}$.

Assuming $Z_{\odot}=0$, the galactocentric radial and azimuthal velocities can
be found from

$\displaystyle v_{R}$ $\displaystyle=$ $\displaystyle-u\cos\phi+v\sin\phi$
$\displaystyle v_{\phi}$ $\displaystyle=$ $\displaystyle u\sin\phi+v\cos\phi$
(12) $\displaystyle v_{z}$ $\displaystyle=$ $\displaystyle w\,,$

where $\phi$ is the galactocentric azimuth, taken as positive in the direction
of galactic rotation:

$\phi=\arctan\left(\frac{Y}{-X}\right)=\arctan\left(\frac{d\cos b\sin
l}{R_{\odot}-d\cos b\cos l}\right)\,,$ (13)

making our $(R,\phi,z)$ galactocentric cylindrical coordinates a left-handed
system.

###  3.3 Propagation of uncertainties

Figure 9: The relative distance uncertainties of the individual sources, as a
function of distance, for the RGB sample with radial velocities. The blue
curve shows the median in bins of 50 pc.

To estimate the uncertainties in positions and velocities from the formal
errors in astrometry, we choose to concatenate the Jacobian matrices of the
consecutive transformations necessary to move from the initial reference frame
to the desired one. In doing so, we are implicitly linearizing the functions
that allow us to convert astrometry into positions and velocities. In other
words, we are simplifying the sequence of non-linear transformations (e.g.,
see eqs. 7,3.2, and 13) that convert the coordinates in one reference frame,
$\@vec{x}_{1}$, to the coordinates in another frame, $\@vec{x}_{2}$, with a
single matrix product of the form

$\@vec{x}_{2}=J\@vec{x}_{1},$ (14)

by taking only the linear term of the Taylor expansion of the transforming
functions. In this case, $J$ corresponds to the Jacobian matrix of the
functions used to transform from one frame to the other:

$J_{i,j}=\frac{\partial f_{i}(\@vec{x})}{\partial x_{j}},$ (15)

where $f_{i}$ is the function that calculates the $i^{\mathrm{t}h}$ component
of the vector $\@vec{x}_{2}$ in the desired coordinate frame given the vector
$\@vec{x}_{1}$ in the original coordinate frame.

In practice, we do the following: first, we construct the covariance matrix in
the frame of the Gaia astrometry using the proper motion and line-of-sight
velocity uncertainties, and their associated correlations. As said above, we
neglect the uncertainties in sky position and consequently all the terms
involved are set to zero. At the same time, we also set to zero the
uncertainty in parallax and its correlations, replacing it with the
uncertainty in our distance estimate, $\sigma_{d}$ (discussed below). The
resulting initial covariance matrix for any source is, thus,

$\Sigma_{ICRS}=\begin{pmatrix}0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&\sigma_{d}^{2}&0&0&0\\\
0&0&0&\sigma_{\mu_{\alpha}^{*}}^{2}&\rho_{\mu_{\alpha}^{*}\mu_{\delta}}\sigma_{\mu_{\alpha}^{*}}\sigma_{\mu_{\delta}}&0\\\
0&0&0&\rho_{\mu_{\alpha}^{*}\mu_{\delta}}\sigma_{\mu_{\alpha}^{*}}\sigma_{\mu_{\delta}}&\sigma_{\mu_{\delta}}^{2}&0\\\
0&0&0&0&0&\sigma_{v_{los}}^{2}\\\ \end{pmatrix},$ (16)

where $\rho_{\mu_{\alpha}^{*}\mu_{\delta}}$ is the correlation coefficient
between the proper motion components, as given in the Gaia source table.

Then, with the Jacobian calculated as in Eq. 15, we obtain the covariance
matrix in heliocentric Cartesian coordinates (eq. 6 and eq. 7):

$\Sigma_{2}=J\Sigma_{1}J^{T}\,,$ (17)

After that, we use again equations 17 and 15 to move, successively, first into
galactocentric Cartesian (eq. 6 and eq. 10) and, finally, to galactocentric
cylindrical (eqs. 3.2 and 13). At each step, we obtain a full 6x6 covariance
matrix that encodes not only the estimated uncertainty in each quantity along
its diagonal, but also their correlations. The transformation itself
contributes significantly to some of these correlations like, for instance,
the correlation between the $u$ and $v$ components of the velocity that arises
naturally from using the distance in eq. 7. Nonetheless, the correlation
between Gaia measurements are also a source of correlations regardless of the
coordinate frames. However, it is the correlations introduced by the
coordinate transformations, in concert with the distance uncertainties, that
dominate the final correlations between the velocity components, and these are
highly direction dependent. (See also the following section for further
discussion.)

Alternatively, we could have chosen to estimate the uncertainties by randomly
sampling new ”fake” observables based on the covariance matrix in the initial
Gaia frame. However, there are a couple shortcomings in this approach: i) it
requires a large amount of samples to obtain a robust estimation of the
uncertainties, which can be computationally expensive, ii) we do not know the
true error distribution function since the formal errors are estimated from
the observables themselves, and iii) we do not have access to the full
posterior distribution function for the distance estimators. While point (i)
can be dealt with some patience, points (ii) and (iii) force us to make
similar assumptions to those made for the strategy described above, thus
limiting the usefulness of this alternative approach.

As a consequence of our decision to use the Jacobians to propagate the
uncertainties, we are implicitly assuming that the uncertainties of the
measured (input) quantities are symmetric. For the proper motions and line-of-
sight velocities this is satisfied, their errors being well described by
Gaussian distributions (Gaia Collaboration et al., 2022b). However, as
described in CBJ2021, the probability distribution functions of the distances
are typically skewed, and this is also true of photometric distances, such as
those used for the Cepheids. For the CBJ2021 distances, we render the distance
uncertainties symmetric by taking the mean of the distances of the provided
16th and 84th percentiles from the median distances, that is, we take as the
individual distance uncertainty:

$\sigma_{d}=(\mathtt{r\\_hi\\_photogeo}-\mathtt{r\\_lo\\_photogeo})/2$ (18)

where $\mathtt{r\\_hi\\_photogeo}$ and $\mathtt{r\\_lo\\_photogeo}$ are
respectively the 84th and 16th percentiles of the “”photo-geometric” distances
provided by CBJ2021. Figure 9 shows how the distance uncertainties in our
sample varies with distance for the subset of RGB stars with radial
velocities. We note that the median relative uncertainties are below 20%. The
distance uncertainties for the OB sample with line-of-sight velocities follows
the same trend, though covering a much smaller range of distances.

Using this distance estimator comes at the cost of loosing the correlations
between the distances and the proper motions. Such correlations are not
expected in purely photometric based distances, as is the case for the
Cepheids, but will be present in as much as the distance is informed by the
parallax for the other distance estimates. For the CBJ2021 distances, used for
the OB and RGB samples, it will be the relatively nearby sources, which
generally have small relative distance uncertainties, that will be constrained
by the parallaxes, while the more distant sources with larger uncertainties
will be primarily constrained by the photometry and the CMD prior, so will be
only weakly correlated with the proper motions. Therefore, to propagate the
uncertainties on the velocities, we only concern ourselves with the
correlations between the proper motion components in the initial covariance
matrix (eq. 16).

Figure 10: The median radial (blue) and tangential (orange) velocity
uncertainties of the individual sources, as a function of distances. The
shaded areas show the range between the 16th and 84th percentiles. The upper
plot is for the OB stars, while the lower plot is for the RGB stars, for the
samples with full velocity information.

We note that, since the tangential (perpendicular to the line-of-sight)
velocity components are dependent on both distance and proper motion, their
uncertainties will be perfectly correlated with the distance uncertainties. As
the line-of-sight velocity does not suffer from correlations with the
astrometry, and may be of a different magnitude with respect to the
uncertainty in the tangential velocity, we should expect that the correlations
and uncertainties of our velocity components will have a strong directional
dependence that can potentially introduce false signals or patterns in our
maps. Figure 10 shows the median and range of the velocity uncertainties for
the OB and RGB samples as a function of distance. We note that the velocity
uncertainty of the OB stars with full velocity information, being limited to
within a few kiloparsecs of the Sun, will be dominated by the uncertainty in
the line-of-sight velocities. On the other hand, the uncertainties for both
velocity components of the RGB stars are quite comparable to about 2 kpc,
beyond which the tangential component then dominates the uncertainty.

Finally, we also note that our distance uncertainties will again introduce
further direction dependencies in the errors and correlations in the final
transformation to velocities in galactocentric cylindrical coordinate, via
equations 3.2 and 13. Indeed, errors in $\phi$ can become quite large near the
galactic centre. Also, by performing the coordinate transformation to
($v_{R},v_{\phi}$) above, we have made ourselves vulnerable to systematic
error in our assumed values of $R_{\odot}$ and ${\bf{v}}_{\odot}$. We explore
the effect of these kinds of uncertainties in the following section.

###  3.4 Effect of systematic errors

Figure 11: Systematic velocity trends induced by the simulated uncertainties
in our mock catalogue. _Upper left panel:_ map of the observed $V_{R}$
assuming a relative parallax uncertainty = 0.2. _Upper right panel:_ same as
left panel, but showing the observed $\Delta V_{\phi}$. _Lower left panel:_
the observed $V_{R}$ assuming that the true value of $R_{\odot}$ is 200 pc
larger than the adopted value, together with the distance uncertainty model
explained in the text. _Lower right panel:_ the observed $V_{R}$ assuming that
the true value of the solar velocity component $V_{\phi,\odot}$ is 10 km/s
larger than the adopted value, together with the distance uncertainty model
explained in the text.

As mentioned above, while the formal errors are quite small there may be
significant systematic errors in our assumed values for $R_{\odot}$ and
${\bf{v}}_{\odot}$. In addition, as we will demonstrate below, even random
distance errors can introduce systematic errors in the mean galactocentric
velocity components. To investigate the possible effects that such errors can
introduce, we construct a mock catalogue from a rather artificial distribution
of stars: We uniformly populate a disk of stars centred on the Sun’s position
with a radius of 8 kpc, add a Gaussian velocity dispersion and an azimuthal
velocity with respect to the galactic centre at $R_{GC}$ that is consistent
with what we observe for the RGB sample (see section 5), but assuming no mean
radial motion (i.e. $v_{R}=0$). We also assume a motion of the observer at the
Sun and from the relative motions derive the proper motion and line-of-sight
velocities, to which we add fractional uncertainties of 0.01 in the proper
motions and 0.1 in the radial velocities. We then re-derive the
($v_{R},v_{\phi}$) velocity components from the observed (noisy) proper
motions, distances and line-of-sight velocities, and construct maps of the
observed mean velocity field, as described in the section 5.

To model the effect of distance uncertainties we add a 20% gaussian
uncertainty to the true parallax, then take the inverse of the observed
parallax as our distant estimate. We note that this is much larger than the
actual uncertainties of our dataset, which only reach a level of 20% relative
uncertainty in the distances at about 10 kpc from the Sun, but is used here
simply for the purpose of illustration. However, like the actual distance
errors, the probability distribution function of the distance from our
inverse-parallax distance estimate is skewed toward larger distances. More
importantly, for our rather artificial uniform disc of mock stars, the mean
estimated distance will be slightly systematically underestimated for
distances less than 8kpc, but strongly overestimated for distances beyond 8
kpc. These mean biases in the distances of our sample introduce systematic
motions in the inferred mean velocity field ($V_{R},V_{\phi}$), as shown in
the upper two panels of figure 11. In the following, we will refer to those
biases as ’systematic errors’, in the sense that the simulated uncertainties
in our experiment systematically induce artefacts in the observed trends, as
shown below.

Assuming a more realistic model for our parallax uncertainties of
$\sigma_{\varpi}/\varpi=0.02d_{\mathrm{k}pc}$ (i.e.
$\sigma_{\varpi}=20\mu{\mathrm{a}s}$), we investigate the effect of assuming
erroneous values for $R_{\odot}$ and the velocity of the Sun with respect the
values used to generate the mock observations. We first introduce a systematic
error in our assumed value for $R_{\odot}$, taken to be 200 pc larger than the
distance to the Galactic centre used to assign the mean rotational velocities
to the mock stars, when we transform the observed proper motions, distances
and line-of-sight velocities to ($v_{R},v_{\phi}$) velocity components. The
resulting inferred $V_{R}$ velocity field is shown in the lower left panel of
figure 11. Similarly, assuming a correct value for $R_{\odot}$ consistent with
the mock velocities, but assuming an incorrect velocity for the Sun’s
$v_{\phi}$ component, we see a similar but inverted pattern in the inferred
$V_{R}$ velocities, with the Sun lying on an axis of symmetry of the inferred
$V_{R}$ velocity field. (Lower right panel of figure 11.)

Another important feature in the velocity maps to be noted are that the
velocities beyond the actual distance limit of our mock sample (8 kpc) are
very strongly biased. This is an unavoidable feature of any magnitude limited
sample, where the real density of the sources will sharply decrease beyond
some limiting distance. Beyond this distance the sources used to estimate
velocities will have systematically over-estimated distances, and consequently
overestimated velocities (in modulus) perpendicular to the line-of-sight. This
limiting distance, beyond which the inferred velocities cannot be trusted, can
vary significantly for different lines-of-sight, thanks to the effects of
interstellar extinction. In any case, one should resist giving astrophysical
significance to features at the edge of velocity maps.

While our assumed systematic errors in the above discussion may be larger than
we expect, the purpose of this discussion is to warn anyone interpreting
features or patterns seen in velocity maps in galactocentric cylindrical
coordinates. As in coordinate space, one should be most suspicious of any
patterns in the kinematics that show any symmetry with respect to the Sun’s
position.

##  4 Coordinate maps

Figure 12: Spatial distribution in the galactic plane, in heliocentric
coordinates, of the selected OB stars (left panel) and RGB stars with line-of-
sight velocities (right panel). For the OB stars, we perform a bivariate
kernel density estimation using an Epanechnikov kernel with a smoothing length
of 200 pc. For the RGB stars, we construct a 2-dimensional histogram, using
cells of 100x100 pc. The galactic centre is toward the right, shown with a
cross in the right hand plot, with galactic rotation being clockwise. Sun’s
position is shown with a black dot. A logarithmic stretch is used for the map
of the RGB stars to enhance lower density regions.

In the present Section, we map the spatial distribution of the OB and RGB
stars, open clusters, and Cepheids samples described in Section 2.
Additionally, we present a comparison with some models available in the
literature, and discuss the observed similarities and differences.

The left panel of Figure 12 shows the distribution of the OB stars in the XY-
plane of the Galaxy. Far from being homogeneous, the distribution of the OB
stars is highly structured, and has numerous regions where the stellar density
is markedly higher than in other ones. Such high-density regions are not
randomly distributed, but appear to be organised in spiral arm segments.
Specifically, we can identify three quasi-diagonal segments crossing the left
panel of Figure 12: from left to right, we discern the Perseus arm, the Local
(Orion) arm, and an inner stripe corresponding to the Sagittarius-Carina and
(possibly) the Scutum arms. Figure 12 (left panel) maps with unprecedented
detail the structural features of the OB stellar population, especially within
about 3 kpc from the Sun. Beyond this distance the distribution becomes
increasingly dominated by radial features produced by foreground extinction.

The right panel of Fig. 12 shows the spatial distribution of the 5.7M RGB
stars with line-of-sight velocities in the Galactic plane. Again, as for the
OB sample, radial ”shadow cones” from foreground extinction are clearly
visible in the RGB sample, though with higher angular frequencies, as the
density here is not derived using a smoothing kernel. Aside from this
difference, and in contrast to the OB stars, the RGB sample exhibits a smooth
spatial distribution, as expected from a dynamically old stellar population.
No clear spiral structure is apparent from the stellar counts, possibly due to
the fact that the giant sample contains typically old stars (compared to the
other populations considered in this work). However, we note that a density
enhancement is present toward the Galactic centre, presumably due to the
Galactic bar. Additionally, we note that an overdense region is apparent at
$x\simeq$ 2-3 kpc, running across a range of y. This is due to the combination
of two main factors: (i) we are using a magnitude limited sample, implying
that the density decreases with heliocentric distance; (ii) the intrinsic
distribution of RGB stars in the Galactic disc is increasing toward the inner
parts of the Galaxy (as expected from an exponential disk).

Figure 13: Overdensity map of the OB stars over plotted with the positions of
the open clusters younger than 63 Myr and with $n_{0}>5$, plotted with filled
circles whose size is proportional to $\sqrt{n_{0}}$. The cross indicates the
position of the Sun.

To better explore the Galactic spiral structure, we apply to our OB sample the
same approach adopted by Poggio et al. (2021) to map the stellar over-density,
defined as

$\Delta_{\Sigma}(X,Y)=\frac{\Sigma(X,Y)\,-\langle\,\Sigma(X,Y)\,\rangle}{\langle\,\Sigma(X,Y)\,\rangle}\,$

where the local surface density $\Sigma(X,Y)$ and the mean surface density
$\langle\,\Sigma(X,Y)\,\rangle$ are constructed using kernel density
estimators with bandwidths of 0.3 and 2 kpc, respectively, adopting an
Epanechnikov kernel. The resulting map is shown in Figure 13. As we can see,
the red diagonal stripes in Figure 13 correspond to the segments of the
nearest spiral arms, consistently with the features identified in the left
panel of Figure 12. Based on Figures 12 and 13, the emerging picture of the
Galactic spiral structure is in good agreement with the one found in Poggio et
al. (2021) (see their Fig. 1B and 1C).

The same Figure 13 shows a comparison between the OB over-density map and the
young ($<63$ Myr) and bright (i.e. with at least 5 members brighter than
$M_{G}=0$) open clusters (OCs). Distances to the OCs were obtained by
inverting the median parallax of each cluster, as described above in section
2.1. As we can see in Figure 13, there is a good agreement between the open
clusters distribution and the spiral structure mapped by the OB sample.

Figure 14: Heliocentric coordinates of the clusters younger than $\log t$=7.6
(63 Myr). Thick grey lines are the spiral arm model of Reid et al. (2019), and
the dashed line is the trace of the Perseus arm modelled by Levine et al.
(2006). The bars represent the 1-$\sigma$ uncertainty on the distance, taking
into account statistical and systematic parallax errors.

Figure 14 shows the spatial distribution of the clusters younger than $\log
t$=7.8 ($\sim$63 Myr). Although the young clusters clearly trace multiple
elongated structures, evocative of arms and inter-arm regions, their
distribution is not continuous. The young OCs alone do not constitute a
sufficient sample to clearly define the main spiral arms. The most striking
difference when comparing this distribution to the spiral arm model of Reid et
al. (2019) (shaded grey arms in Fig. 14) is that the Perseus arm appears
interrupted for two kiloparsecs. This discontinuity has been observed before
in cluster distributions (e.g. Cantat-Gaudin et al., 2020) or CO clouds (Peek
et al., 2022). On the other hand, the model of Levine et al. (2006) traces
very nicely the orientation of the Perseus arm in the upper main sequence
stars (see also Poggio et al., 2021), and appears to be in reasonable
agreement with the distribution of the OCs as well (see the orange dashed line
in Fig. 14). According to this model, the two groups that Reid et al. (2019)
consider to define a low-pitch-angle Perseus arm (one in the second and one in
the third Galactic quadrant) would in fact belong to two different arms.

Figure 15: Wavelet transformation of the Cepheids with age ¡ 200 Myr. Black
dots in the left panel show the positions of the single sources, while the
right panel shows, on a larger scale, an overlay of the model from Taylor &
Cordes (1993) (solid lines) and the one from Levine et al. (2006) (dashed
line).

Finally, we map the spatial distribution of Cepheids younger than 200 Myr
using a wavelet transformation (WT). The left panel of Figure 15 shows a
comparison between the single sources (shown as black dots) and the WT
coefficients (colored map, see technical details in Ramos et al., 2018; Poggio
et al., 2021). The right panel of Figure 15 shows a comparison between the
Cepheids WT and some models available in the literature. Solid lines show the
4-armed model of Taylor and Cordes (Taylor & Cordes, 1993); in the region
where the Cepheids are more abundant ($x<5$ kpc), the Sagittarius Carina arm
appears to be in good agreement below approximately $y\simeq$2.5 kpc. At
approximately $x\simeq 2.5$ kpc and $y\simeq$ 5 kpc, the overdensity of
Cepheids seem to diverge with respect to the Taylor and Cordes model of the
Sag-Car arm in this direction. Meanwhile, in the outer regions of the Galaxy,
the orientation of the Perseus arm in the Cepheids seems more consistent with
the Levine model (Levine et al., 2006) than either the Taylor and Cordes or
Reid models. This outer arm, based on HI data, is remarkably traced by the
Cepheids out to a galactocentric radius of at least 16 kpc. Unfortunately the
Cepheids are too sparse to trace the weaker Local (Orion) arm.

##  5 Velocity maps

This section studies the kinematics of RGB and OB stars in the disc of the
Milky Way to highlight the effects of the disc asymmetries on velocities. We
use the selections of OB and RGB stars within $|z|\leq 0.3$ (77 659 stars) and
$|z|\leq 1$ kpc (5 730 578 stars) of the Galactic plane, respectively. We
first describe the construction of the maps of the mean velocity and velocity
dispersion from the individual velocities derived in section 3, and analyse
the resulting velocity fields of the RGB sample, which cover a much larger
extent of the disc than the OB sample, discussed at the end of this section.
Table 4 summarising our notation is provided for convenience.

Table 4: Notation and Nomenclature variables |
---|---
$\varpi,\@vec{\mu}$ | astrometry
$\sigma_{\varpi},\sigma_{\mu_{\alpha}},\sigma_{\mu_{\delta}},$ | their uncertainties
$v_{\mathrm{los}}$ | line-of-sight (los) velocities
$d,\sigma_{d}$ | distance and related uncertainty
$x,y,z$ | heliocentric Cartesian coordinates
$X,Y,Z$ | galactocentric Cartesian coordinates
$u,v,w$ |  $x,y,z$ velocity components
$R,\phi,z$ | LH galactocentric cylindrical coordinates
$v_{R},v_{\phi},v_{z}$ |  $R,\phi,z$ velocity components
$\sigma_{v_{R}},\sigma_{v_{\phi}},\sigma_{v_{z}}$ | their uncertainties
$V_{R},V_{\phi},V_{z}$ | mean $v_{R},v_{\phi},v_{z}$ velocities
$\sigma_{V_{R}},\sigma_{V_{\phi}},\sigma_{V_{z}}$ | their uncertainties
$V_{R}$ | radial velocity
$V_{\phi}$ | azimuthal velocity
$\overline{V}_{\phi}(R)$ | mean azimuthal velocity at $R$
$\sigma^{*}_{R},\sigma^{*}_{\phi},\sigma^{*}_{z}$ | velocity dispersions in $R,\phi,z$ directions

###  5.1 Construction and analysis of velocity maps

We built maps of the ordered and random motions for the radial, azimuthal, and
vertical velocity components. We designed these maps with a constant $100$ pc
resolution, with $341\times 341$ pixels. These characteristics result from
empirical choices to have sufficient numbers of stars per bin for the current
analysis. To perform robust derivations of the velocity per pixel, we
considered only the cells with a minimum of $20$ stars, and masked all others.
This gives grids with a median and a maximum numbers of $58$ and $267$ stars
per pixels for OB stars, respectively, and $152$ and $2\,541$ stars per pixel
for RGB stars, respectively.

At a given heliocentric position $(x,y)$ corresponds a cell $j$ of $100\times
100$ pc in our three velocity components $k=R,\phi,$ or $z$ maps, which
contains $N_{*}$ stars. We estimate the stellar mean velocity $V_{k}$ and its
associated dispersion101010Most inconveniently, $\sigma$ is the traditional
notation for velocity dispersion, but also that for Gaussian uncertainties. To
avoid potential confusion, we add a superscript $*$ to the stellar velocity
dispersion. $\sigma^{*}_{k}$, by optimising the log-likelihood of the
distribution of “observed” velocities $v_{k,j}$ (with uncertainties
${\sigma_{v}}_{k,j}$) of the stars located within the $j$-cell. We assumed
Gaussian uncertainties on the individual independent $v_{k,j}$, so that the
negative log-likelihood to minimise is:

$\mathcal{L}(V_{k},\sigma^{*}_{k})=\frac{1}{2}\sum^{N_{*}}_{j}\left(\ln({\sigma^{*}_{k}}^{2}+{\sigma_{v}}_{k,j}^{2})+\frac{(v_{k,j}-V_{k})^{2}}{{\sigma^{*}_{k}}^{2}+{\sigma_{v}}_{k,j}^{2}}\right).$
(19)

Appendix B.1 details the derivations of the uncertainty ${\sigma_{V}}_{k,j}$
of our mean velocities ${V}_{k,j}$, which also demonstrates that under certain
conditions we can approximate them as
$\sigma_{V_{k}}=\sigma^{*}_{k}/\sqrt{N_{*}}$.

Figure 16: Velocity maps of RGB stars. Left and right panels are for the
ordered and random motions, respectively. From top to bottom, panels show the
radial, azimuthal, and vertical velocity components. The origin, emphasised by
the black dot, indicates the position of the Sun, and a plus symbol indicates
the Galactic Centre. Dashed circles represent constant distance to the
Galactic centre at $R=5$ and $R=10$ kpc. The velocity ranges have been chosen
to enhance contrasts to help identify by eye regions where streaming motions
occur. Left panels of Fig. 24 show the associated uncertainty maps. We detail
the construction of these maps in Sect. 5.1. These velocity and velocity
dispersion maps can be found in CDS in FITS format.

Figure 16 shows the resulting three component velocity fields for the sample
of RGB stars (left panels), and the velocity dispersions (right panels). We
optimised the displayed velocity ranges to help identify by eye regions where
streaming motions occur. Left panels of Fig. 24 shows the associated
uncertainty maps for the velocities.

The $V_{R}$ map shows a remarkable bisymmetric feature on both sides of the
GC, with negative and positive values on each side of the apparent bar’s major
axis. This quadrupole feature is a characteristic of the mean inward motion
down to $\sim-40$ $\rm\,km\,s^{-1}$ ($y>0$) and the mean outward motion up to
$\sim 45$ $\rm\,km\,s^{-1}$ generated by the Galactic bar. Gaia DR3 allows us
to confirm the bar quadrupole $V_{R}$ pattern identified in Bovy et al. (2019)
and Queiroz et al. (2021), who used thousands of stars, but using Gaia DR2 and
EDR3 astrometry (Gaia Collaboration et al., 2018a, 2021b) together with APOGEE
line-of-sight velocities (Majewski et al., 2017; Abolfathi et al., 2018).
Section 5.3 and Appendix C describe in detail this quadrupole pattern. We also
note that the RGB sample contains enough stars to apparently provide us with
mean velocity estimates beyond the Galactic Centre (GC), though these should
be interpreted with caution.

The $\sigma^{*}_{R}$ map shows a bisymmetric pattern as well, but different
from $V_{R}$, with larger amplitudes which are aligned with the direction
where $V_{R}$ changes its sign along the bar major axis, and lower along a
perpendicular direction. Here again, the GC is the node of the quadrupole
feature. In addition to the central quadrupole, streaming motions in $V_{R}$
also occur at larger radii, e.g. $R=6.5$ kpc ($x\sim 1.5$ kpc), which shows
$V_{R}$ being larger for $|y|<2$ kpc than at smaller and larger azimuths. At
$R\sim 10$ kpc ($x,y\sim-1.5,+1$ kpc), $V_{R}$ shows a clear change of sign
with respect to azimuth.

The distribution of the azimuthal velocity $V_{\phi}$ is elongated in the bar:
within the central 5 kpc, the rotation at a given radius is slower along the
apparent bar axis than perpendicular to the bar axis. The $\sigma^{*}_{\phi}$
map also seems to exhibit a bisymmetry in the bar region, but rotated by about
$45\degr$ with respect to that seen in $\sigma^{*}_{R}$. The azimuthal random
motion appears smaller when the radial random component is larger. In other
words, the planar velocity dispersion is highly anisotropic in the bar region.

The vertical velocity $V_{z}$ is mostly positive. Unlike the other two
velocity components, it does not show any pattern linked to the presence of
the Galactic bar. A streaming is observed around $R=7$ kpc, as $V_{z}$ shows
among the lowest values on one side ($1\leq x\leq 4,y<-2$ kpc), while being
positive on the opposite side with respect to $y=0$. However, this pattern,
being symmetric about the x-axis, may be an artefact of systematic errors.
Vertical motions are also larger with galactocentric radius toward the
galactic anticentre, clearly showing the kinematic signature of the warp of
the Galactic disc beyond, in agreement with previous results (Gaia
Collaboration et al., 2018b; Romero-Gómez et al., 2019; Poggio et al., 2018,
2020; López-Corredoira & Sylos Labini, 2019). We do not observe any
bisymmetric feature in the vertical dispersion, which is larger along $l\sim
0$. It should be noted that the vertical velocity component is more sensitive
to systematic errors as, our sources being near the galactic ($b=0$) plane,
this component is predominantly manifested in the tangential velocities that
are sensitive to distance errors (see Fig. 10).

In summary, these maps detect the significant signature of the Galactic bar of
the MW in the inner disc. We also find some streaming signatures of $V_{R}$ at
larger galactocentric distances, which might be associated with corotation or
the outer Lindblad resonance. Section 5.4 discusses the case of the young
stellar populations and the similar maps for the OB star sample.

###  5.2 Analysis of radial profiles

Figure 17: Velocity profiles of the Milky Way RGB and OB stars (black and blue
lines, respectively). Top panel: Rotation curves. Shaded area represent the
uncertainties. Bottom panel: Radial, azimuthal and vertical velocity
dispersions (solid, dashed and dotted lines, respectively). The data shown in
these figures can be found in CDS.

From these maps, we inferred the average axisymmetric variation of velocities
following the same procedure as in Gaia Collaboration et al. (2021c), where
for each bin in radius of 200 pc size, the median value of the cells in the
radial bin define the azimuthally-averaged velocity $\overline{V}_{\phi}(R)$
at that radius. We discarded the radial bins where the number of pixels from
the maps is less than 5 pixels. Bootstrap resamplings were performed at each
radius to define the velocity uncertainties, measured at the 16th and 84th
percentiles of the velocity distributions. In Fig. 17, we present the Galactic
rotation curves, as well as the velocity dispersion profiles for the RGB and
OB samples. Despite the asymmetries observed in the velocity field, especially
at radii $R<5$ kpc, the rotation curve of the RGB stars is regular. It
smoothly increases like a solid body up to $\overline{V}_{\phi}\sim 220$
$\rm\,km\,s^{-1}$ to $R\sim 6$ kpc, then remains constant out to the last
measured radius. The amplitude of the curve for the younger OB stars is larger
than for RGB giants, because of their smaller asymmetric drift, and is
decreasing with galactocentric radius, rotating on average 17
$\rm\,km\,s^{-1}$ faster than the RGB stars over the radial range $R=6.5-10$
kpc. These average rotation curves have been subtracted from the $V_{\phi}$
maps to produce residual velocity fields
$V_{\phi}(x,y)-\overline{V}_{\phi}(R)$, useful for evidencing velocity
streaming in the outer disc (Sect. 5.4).

Figure 18 shows a comparison between the azimuthal velocity of the OB stars,
the open clusters, the Cepheids, and the RGB giants with respect to
galactocentric radii. The median velocities for Cepheids and open clusters are
calculated using overlapping radial bins of 1 kpc step and 2 kpc width. As
expected, there is a good agreement between the OB stars and the OCs. The
Cepheids exhibit a mean azimuthal velocity similar, but slightly slower, to
the OB and the open clusters, but reach significantly larger radii. Finally,
the mean azimuthal velocity for the RGB giants is systematically lower than
the other tracers, as expected for asymmetric drift for an older population
with a larger velocity dispersion.

The radial and azimuthal velocity dispersion profiles of RGB stars show two
distinct regions. Within the inner part of the bar ($R\sim 2.5$ kpc), the
profiles are shallow and the dispersions comparable ($\sim 75-80$
$\rm\,km\,s^{-1}$), though we note that the $(x,y)$ map of the dispersion
(Fig. 16) in this region shows very strong asymmetries. Beyond 2.5 kpc,
$\overline{\sigma}^{*}_{z}<\overline{\sigma}^{*}_{\phi}<\overline{\sigma}^{*}_{R}$
and the profiles continuously decrease for all three components. The radial
profiles of the velocity dispersion $\overline{\sigma}^{*}_{k}$ of the OB
stars do not vary strongly with radius, showing average radial, azimuthal and
vertical dispersions of $14.2,9.4$ and $6.2$ $\rm\,km\,s^{-1}$, respectively,
in very good agreement with the observed random motions of gas within $R=8$
kpc ($4-9$ $\rm\,km\,s^{-1}$, Marasco et al., 2017). The axis ratios of the
velocity ellipsoid as averaged from these profiles for $R>3$ kpc are
($\overline{\sigma}^{*}_{\phi}/\overline{\sigma}^{*}_{R},\overline{\sigma}^{*}_{z}/\overline{\sigma}^{*}_{R},\overline{\sigma}^{*}_{z}/\overline{\sigma}^{*}_{\phi})=(0.81,0.60,0.75)$
for the RGB stars, and $(0.66,0.44,0.66)$ for OB stars.

Figure 18: Comparison of rotation curves of the different samples. Mean
azimuthal velocity for the open clusters (dark blue) and Cepheids (yellow) as
a function of galactocentric radius (calculated as explained in the text),
compared to the profiles of the OB (light blue) and RGB stars (turquoise).
Coloured points show the single sources for the Cepheids (black dots) and OCs
(red points).

###  5.3 The kinematics of the bar

In this Section, we estimate some fundamental parameters of the Galactic bar,
guided by results obtained from Section 5.1 and with the help of a mock galaxy
from a numerical test particle simulation of a barred galaxy. We want to make
clear that this is not a made-to-measure, customised simulation that intends
to quantitatively reproduce the observed dataset of RGB stars. The simulation
has to be seen as a simple diagnostic tool for qualitative comparison to the
observations, and the results presented below as tentative possibilities. We
refer the reader to App. C) for the description and analysis of the
simulation. We thus applied the same recipes as those described in App. C) to
find, by analogy, the orientation and pattern speed of the Galactic bar, as
well as the location of the outer Lindblad resonance.

We start to fit the non-uniformity of the kinematics in the bar region of the
RGB sample, assuming that the Galactic bar perturbs velocities by adding a
bisymmetric component to the axisymmetric motions. This should apply to most
of the ordered and random motions, because all of them but $v_{z}$ were shown
to be similarly structured in the bar region (Sect. 5.1). We therefore
performed a simple Fourier decomposition up to second order to characterise
the axisymmetric and the bisymmetric components seen in the maps. This
approximation of $V_{R}$ and $V_{\phi}$ is referred to as $V_{R,\rm mod}$ and
$V_{\phi,\rm mod}$, and is given by:

$V_{R,\rm
mod}(R,\phi)=\overline{V}_{R}(R)+A_{R}(R)\cos\left(2(\phi-\phi_{R}(R))\right),$
(20)

where $\overline{V}_{R}(R)$ is the axisymmetric mean value, and $A_{R}$ and
$\phi_{R}$ are the amplitude and phase angle of the bisymmetric Fourier
harmonics of $V_{R}$ respectively. Another parameter of interest is the
scatter in the modelling, $V_{R,s}$, which absorbs all other asymmetric
departures from the bisymmetry in the velocity fields. We use an analogous
expression for $V_{\phi,\rm mod}$.

Figure 19: Results of the bisymmetric model of the $V_{R}$ and $V_{\phi}$
velocity maps. Top: Amplitudes $A_{R}$ and $A_{\phi}$ of the bisymmetry. The
axisymmetric component $\overline{V}_{R}$ is also shown. Bottom: Phase angle
of the perturbation with respect to the Sun-GC direction. Shaded area
represent the uncertainties. The horizontal dotted line marks the estimated
bar orientation angle. The vertical dashed line and purple area show the
estimated location of the bar corotation.

We performed Bayesian inferences of the model through Markov Chain Monte Carlo
(MCMC) fits, using the Python library emcee (Foreman-Mackey et al., 2013). As
this model best applies to regions where the kinematics is perturbed by a
second order perturbation, we restricted the analysis to $0\leq R\leq 7$ kpc,
and considered a radial bin width of $200$ pc. Defining the residual velocity
as $V_{R,\rm res}=V_{R}-V_{R,\rm mod}$, the conditional likelihood function
for each radial bin is expressed by:

$\mathcal{L}(\overline{V}_{R},A_{R},\phi_{R},V_{R,s})=-\frac{1}{2}\left(n_{\rm
pix}\ln(2\pi)+\sum^{n_{\rm pix}}\left(V_{R,\rm
res}^{2}/\xi^{2}+\ln(\xi^{2})\right)\right),$ (21)

and similarly for $\phi$, where $\xi^{2}=\sigma^{2}_{V_{R}}+V^{2}_{R,s}$, and
$n_{\rm pix}$ is the number of pixels inside the corresponding radial bin,
imposing $n_{\rm pix}>25$ as a condition. Radial bins not satisfying this
condition are not considered. We set a number of $32$ walkers and $2000$ steps
in the MCMC fits, which is enough to converge towards robust and stable
solutions. The prior distributions are uniform and span $[-50,50]$
$\rm\,km\,s^{-1}$ for $\overline{V}_{R}$, $[0,300]$ $\rm\,km\,s^{-1}$ for
$\overline{V}_{\phi}$, $[0,50]$ $\rm\,km\,s^{-1}$ for $A_{R}$, $[0,30]$
$\rm\,km\,s^{-1}$ for $V_{R,s}$, and $[-80,110]\degr$ for $\phi_{R}$.
Following prescriptions from App. C, $\phi_{R}$ and $\phi_{\phi}$ are linked
to the direction of the bisymmetric perturbation of density with respect to
the Sun-GC direction in the bar region, and that we call $\phi_{\textrm{b}}$,
by $\phi_{\textrm{b},R}=\phi_{R}-\pi/4$ in the case of the $V_{R}$ model, and
$\phi_{\textrm{b},\phi}=\phi_{\phi}-\pi/2$ in the case of $V_{\phi}$. We quote
the uncertainties on the parameters at the 16th and 84th percentiles of the
posterior distributions.

Figure 19 shows the resulting fits, where it can be seen that the bisymmetry
is strongest at $R=2.5$ kpc in both the azimuthal and radial velocity fields.
While the amplitude of the perturbation decreases beyond 2.5 kpc for $A_{R}$,
it admits other maxima at $R\sim 0.8$ and $6$ kpc for $A_{\phi}$.
Unsurprisingly, this simple model shows that $V_{R}$ is anything but
axisymmetric, as $A_{R}$ exceeds $\overline{V}_{R}$. The axisymmetric radial
velocity is mostly positive, showing a bulk outward motion of 3.6
$\rm\,km\,s^{-1}$ for $R<3$ kpc, and null beyond 3 kpc, on average. As the
radial velocity should be null on average, for a disc nearly relaxed, we
attribute this non-negligible $\overline{V}_{R}$ to the incomplete coverage of
azimuthal angles. In other words, this bulk motion is only representative of
the observed portion of the Galactic disc. Another feature of interest is the
low velocity scatter parameter $V_{R,s}$ and $V_{\phi,s}$, which is $\sim 6$
$\rm\,km\,s^{-1}$ at small radius and decreases to $<1$ $\rm\,km\,s^{-1}$ out
to $R=7$ kpc (not shown in Fig. 19). As for $\overline{V}_{\phi}$, we find it
is very similar to the median rotation curve derived in the previous section.

The orientation of the bisymmetry with respect to the Sun-GC direction is
found to be different in the $V_{R}$ and $V_{\phi}$ models, particularly
beyond $R=3$ kpc as seen in the bottom panel of Fig. 19. At the peak of the
strength at $R=2.5$ kpc, $\phi_{\textrm{b}}$ differ by $\sim 12\degr$.
Interestingly, $\phi_{\textrm{b}}$ in the $V_{\phi}$ model remarkably shows
the same trend as the one seen in the numerical simulation (Fig. 26 of App.
C), where it is seen that it faithfully traces the true bar orientation, even
in the presence of uncertainties. As shown in Fig. 27, the $\phi_{\textrm{b}}$
in the $V_{R}$ model is more affected by the uncertainties. We can thus infer
that the bar angle with respect to the Sun-GC direction is very close to the
value for which $\phi_{\textrm{b},\phi}$ remains flat before the abrupt drop,
that is $-19.2\degr\pm 1.5\degr$ for $1.5<R\leq 2.8$ kpc. Then, as the
location of the minimum of $\phi_{\rm bisym}$ of the $V_{\phi}$ model beyond
the drop of phase and before it rises at larger radii corresponds to the
corotation of the galactic bar in the simulation (see right panel of Fig. 26),
we find by analogy that the range $R=5.2-6$ kpc is hosting the corotation
radius of the Galactic bar, $R_{\rm CR}$. We adopt a conservative value
$R_{\rm CR}=5.4\pm 0.5$ kpc. This is surprisingly also the location where
$\phi_{\rm bisym}$ starts to decrease for $V_{R}$. This coincidence was not
seen in the numerical simulation, maybe because of the lack of spiral
structure beyond the bar in the mock data.

Figure 20: Lindblad resonances of the Galactic bar. The black solid (dotted)
line shows the angular frequency (angular frequency plus half of the epicyclic
frequency) for the RGB stars, while the blue solid line shows the angular
frequency for the OB stars. The horizontal dashed line represents the bar
pattern speed, and the two vertical dashed lines at $R=5.4$ and $9.7$ kpc,
with their uncertainty range in purple, are the location ranges of the
corotation and outer Lindblad resonance, respectively.

In Fig. 20, we show the angular velocity curve,
$\Omega=\overline{V}_{\phi}/R$, and its combination with the epicyclic
frequency for the RGB stars (black) and OB stars (blue, only the angular
velocity). We emphasise here that the epicyclic frequency is derived by
assuming the epicycle approximation, which might not be perfectly correct in a
radial range far from the Solar neighbourhood. Assuming the corotation radius
range obtained above from the phase of the bisymmetry in $V_{\phi}$, the
intersection between $R_{\rm CR}$ and the $\Omega$ curve provides the bar
pattern speed, which is $\Omega_{\rm bar}=38.1^{+2.6}_{-2.}$ $\rm\,km\,s^{-1}$
kpc-1. We can also provide an estimation of the Outer Lindblad resonance
(OLR), which is the galactocentric radius where the fixed pattern speed
intersects the curve $\Omega+\kappa/2$, that is, the position in the disc
where stars rotate slower than the bar pattern and perform two radial
oscillations for one revolution around the Galactic centre ($\kappa$ is the
epicyclic frequency). In this case, the RGB sample provides an OLR of $9.7\pm
0.5$ kpc, outside the Solar radius. A word of caution has to be written here,
because we estimate $\Omega_{\rm bar}$ and the OLR from the azimuthal angular
velocity of the RGB sample, which is not the same as that from the circular
velocity due to asymmetric drift. We note that the OB stars, being younger and
supposedly less affected by asymmetric drift, have an angular velocity
slightly larger than that of RGB stars. By extrapolation to the inner disc,
our $\Omega_{\rm bar}$ value could thus be biased towards lower values and
should be taken as a lower limit.

###  5.4 Kinematics of the outer disk

Figure 21: Velocity maps of OB stars. Left panels show the inferred radial,
azimuthal, and vertical velocity, from top to bottom. The velocity ranges have
been chosen to enhance contrasts. The right panels show the maps of the radial
velocity (top panel), the azimuthal residual velocity (middle panel) and the
vertical velocity (bottom panel), compared to the overdensity contours of the
OB stars (grey shaded areas), calculated as explained in Section 3.1. In all
panels, the position of the Sun is marked by a filled circle. The maps were
obtained considering only stars within $|z|\leq 0.3$ kpc. The velocity and
maps (left hand column) can be found in CDS in FITS format.

In this Section, we present and discuss the velocity maps for the OB and RGB
sample in the outer disk (i.e. with galactocentric radius $R>5$ kpc), and
compare them to the spatial location of the spiral arms in the Galaxy.

Figure 21 shows the velocity maps of the OB stars, with the
$V_{R},V_{\phi},V_{Z}$ maps on the left panels. On the right panels we show
the same maps (with the exception of the middle panel, where $V_{\phi}$ has
been replaced by the residual map of $V_{\phi}$, calculated as explained in
Sec. 5.1), but compared to the spiral arms found in the overdensity (see Sec.
4), overlayed as gray-shaded contours. We note that here, for the OB stars,
the overdensity in the spiral arms and the streaming motions can be studied
using _the same stellar population_ (although, of course, the two samples are
not exactly the same, given that only the OB stars with line-of-sight
velocities have been used to derive the velocity maps). This gives us
confidence that the maps are self-consistent, and that the comparison between
the spatial spiral arms and the corresponding streaming motions is
appropriate.

The top-right panel of Figure 21 shows an alternating positive-negative
pattern in $V_{R}$, the orientation of which is not aligned with the spiral
arms in density. A prominent feature of stars with positive $V_{R}$ is
apparent at approximately $x\approx 0$ kpc, crossing the map almost vertically
and connecting the upper edge of the Local Arm with the lower side of the Sag-
Car arm. The middle-right panel of Figure 21 shows the residual of $V_{\phi}$
with respect to the mean $\overline{V}_{\phi}(R)$. Based on this map, the
stars located just outside and inside the Local Arm move systematically slower
than the mean $\overline{V_{\phi}}$ (yellow-green regions), while stars lying
on the Local Arm appear to move systematically faster than
$\overline{V_{\phi}}$ (red region). It should be noted, however, that the
alignment of these azimuthal streaming motions with the density contours is
not perfect in the lower part of the map. On the other hand, a striking
alignment between the local arm and a systematically positive vertical
velocities $V_{z}$ is shown in the bottom-right panel of Figure 21, indicating
that the OB stars exhibit both in-plane streaming motions as well as vertical
bending waves, all with a relatively short radial wavelength.

Figure 22:  Velocity maps for the RGB sample. The left panel shows the
inferred radial velocity, $V_{R}$. The right panel shows the map for
$V_{\phi}-\overline{V_{\phi}}$, i.e., the residual azimuthal velocity.
Overlaid on both maps is the 2-arm $NIR$ spiral model from Drimmel (2000),
shown as solid lines. The dashed lines correspond to the minimum interarm
density (see text). The position of the Sun is marked by a black dot.

In contrast with the OB sample, the RGB sample is an older, dynamically
relaxed stellar population extending much further from the Sun, allowing us to
trace the outer disc kinematics over a much larger extent than the OB stars.
In order to better highlight the observed features in the outer disc, we mask
out the region $R<5$ kpc and show the $V_{R}$ map and the $V_{\phi}$ residuals
in Figure 22. To compare the velocity maps of the RGB stars to the spatial
position of the spiral arms in the NIR, we make a comparison with the 2-arm
model from Drimmel (2000), shown as solid lines, while the $\pi/2$ phase
shifted geometry, corresponding to the minimum interarm density, is shown as a
dashed curve. While no clear spiral structure is evident in the RGB spatial
distribution (see Section 3.1), the large scale streaming motions in the
radial velocities suggest some possible correspondence, especially along the
inter-arm dashed curve.

Also noteworthy, the $V_{R}$ map (left panel) shows a large positive $V_{R}$
feature in the third quadrant (between $180^{\circ}<l<270^{\circ}$), which
becomes positive at about $l=170\degr$ between about $R=9$ and 11 kpc. This
feature is approximately aligned with an analogous change in the radial
velocities in the inner bar region (see also the top-left panel of figure 16),
and is similar to a feature seen in simulations with a bar near the Outer
Lindblad resonance. In contrast to the $V_{R}$ map, the $V_{\phi}$ residuals
in the right panel of Figure 22 do not show any clear non-axisymmetric
features in the outer disk. There is, however, a noticeable gradient in the
$V_{\phi}$ residuals, about $y=0$, where the residuals are systematically
lower for larger values of $|y|$. Being symmetric about the $x-$axis, this is
likely a systematic in the velocity due to distance uncertainties, as noted in
sec 3.4 (Figure 11). In any case, we remind the reader that features near the
edge of these maps may well be artefacts from oversampling sources with over-
estimated distances, as also discussed in sec 3.4.

In addition to exploring non-axisymmetry in the disc, with the increase in
spatio-kinematic coverage, we can also explore (a)symmetry about the plane of
the Galaxy. In Figure 23 we present velocity difference maps between the $z>0$
& $z<0$ hemispheres, for the RGB sample (top panels) and OB stars (bottom
panels). In each of the components ($k={R,\phi,Z}$) we make velocity maps for
both the upper and the lower hemispheres, then take their difference (i.e.
$\Delta V_{k}=V_{k,Z>0}-V_{k,Z<0}$). As with the other velocity maps, we have
also overplotted the 4-arm spiral model (black) based on near infrared (NIR)
data from Drimmel (2000). Additionally, we plot the $Perseus$ (red), Sag-Car
(purple), and the Local (cyan) arms from Reid et al. (2019). For the RGB
stars, Figure 23(a) shows that in $\Delta{V_{\phi}}$, within a heliocentric
radius of about 3 kpc, there is no significant difference between the upper
and the lower disc, while beyond $R>11$ kpc the disc is rotating faster below
the disc plane ($z<0$) by up to 10$\rm\,km\,s^{-1}$, in the third quadrant and
in the anticentre direction. A similar feature is seen in the $\Delta{V_{R}}$
map, though here in the opposite sense, with more positive radial motion above
the disk plane. By comparing with Figure 22 (left panel), we see that this
feature corresponds with the positive $V_{R}$ feature in the third quadrant
already noted above. It becomes clear that the positive $V_{R}$ in this part
of the outer disc is almost entirely located above the disc plane. It is
interesting to note that, about the plane of the Galaxy, $\Delta{V_{R}}\sim
0$$\rm\,km\,s^{-1}$ in the rest of the map. Finally, in Figure 23(c), we show
the $\Delta{V_{z}}$ map. This is equivalent to mapping breathing modes
(contraction and expansion with respect to $z=0$) as in Widrow et al. (2014).
The absolute amplitude of the $\Delta{V_{z}}$ is about 2-4 $\rm\,km\,s^{-1}$,
i.e., much lower than in the $\phi$ and $R$ components. Nevertheless, again
the outer disc shows a marked feature, though more toward the anticentre.
Also, we note that at least one of the features in the $\Delta{V_{z}}$ map
roughly aligns with the location of the Perseus arm from Reid et al. (2019).

The $\Delta V$ map for the OB stars cover a much smaller extent of the disk.
Indeed we can only construct difference maps not much beyond 1 kpc from the
disc, shown in the lower panels of Figure 23, with the spiral overdensity
contours overlaid as in 21. In Figure 23(d), we note that that residuals in
$\Delta{V_{\phi}}$, are positive in between the arms. This would suggest that
inside the arms, OB stars rotate faster in the upper disc. In Figure 23(e), we
note that $\Delta{V_{R}}$ is generally negative inside the arms, implying that
the OB stars in the upper disc are moving towards the Galactic centre inside
these arms. In Figure 21, we already noted that for the OB population,
$\Delta{V_{R}}$ is generally negative between the arms; the $\Delta{V_{R}}$
map for the OB stars shows that this is stronger in the upper disc. Similarly,
Figure 21 also showed that the $\Delta{V_{\phi}}$ residuals were also higher
in the inter-arm region, and the $\Delta{V_{\phi}}$ map shows this is more
prominent in the upper disc. This suggests there may be some shearing motion
in the $R$ and $\phi$ directions associated with these arms, or that the
spiral perturbation that aligns with the OB population is stronger in the
upper disc. Finally, for completion, we also present the $\Delta{V_{z}}$ map
for the OB stars, in Figure 23(f), showing that there may be an associated
breathing mode, though here we see a net expansion that may be aligned with
the overdensity of the Local arm. In general, the amplitude of the $\Delta
V_{\rm k}$ for the OB stars is much lower compared to the RGB stars.

Figure 23: Heliocentric maps of $\Delta V_{\rm i}$, i.e., ($V_{Z>0}-V_{Z<0}$),
for the components s, $i={\phi,R,Z}$, in the RGB sample (top), and the OB star
sample (bottom). Also overplotted is the near infrared 2-arm spiral model
(black) from Drimmel (2000), as well as the $Perseus$ (red), and the $local$
(cyan) arms model from Reid et al. (2019). The position of the Sun is
indicated by the black dot at $(x,y)=$(0,0).

##  6 Discussion

In section 4 we have summarised the distribution of our four samples in the
$(x,y)$ plane. As expected, the younger samples indeed trace well-known spiral
arm segments, confirming earlier results using Gaia EDR3 astrometry and young
stellar samples selected using NIR photometry (Zari et al., 2021; Poggio et
al., 2021). Segments of the nearest spiral arms are evident, though we have
made no attempt to derive the local stellar density. Doing so would require an
accurate 3D map of extinction, as well as the selection function of the sample
being used, which depends both on the survey selection function as well as the
criteria used to construct the sample of tracers (Rix et al., 2021). For
instance, for compiling our OB and RGB samples we have used the new
astrophysical parameters in Gaia DR3 , which are only available for a fraction
of the stars. As reported in (Creevey et al., 2022; Fouesneau et al., 2022;
Andrae et al., 2022), about half of the brighter stars are missing
astrophysical parameters, which leads to an artificial depression in the
apparent stellar density within about 1 kpc of the Sun.

In contrast, as already discussed in section 4, extinction is responsible for
sampling bias with respect to angle, leading to the evident radial features
centred on the Sun’s position. The radial features so obviously visible in the
distribution of our sources are only weakly visible in the velocity and
dispersion maps in section 5, and are a consequence of biased distance
estimates. The fact that they are not very evident gives us confidence that we
are justified in assuming that our samples are not kinematically biased, and
that the adopted distances are reliable. Nevertheless, as discussed in section
3.4, beyond some limiting distance the majority of stars in a given volume
element will have over-estimated distances, which will in turn lead to the
mean velocity component tangential to the line-of-sight being over-estimated.
This will cause a direction dependent distance limit to the velocity maps,
beyond which systematic errors will dominate. Such systematic errors will show
symmetries with respect to the Sun-galactic centre line. For example, such a
signature is evident in the map of $V_{\phi}-\overline{V_{\phi}}$ (figure 22),
showing large negative azimuthal velocities for distances beyond about 6 kpc
in the directions of $40^{\circ}<|l|<60^{\circ}$. Keeping in mind the
principles noted at the end of section 3.4, we discuss further below the
kinematic maps presented in section 5.

###  6.1 Kinematics

While the optical passband of Gaia limits its mapping capabilities in
configuration space, it nevertheless samples the kinematics of the disc of the
Milky Way to large distances. Here, as is the usual practice, we assume that
our samples are not kinematically biased. This is a reasonable assumption as
long as each sample either covers a limited range of ages, or can be
considered as kinematically relaxed and of common origin. In particular, we
have not made any effort to make a distinction between the young (thin) disc
stars and the old $\alpha$-rich (thick) disc in our RGB sample (Gaia
Collaboration et al., 2022c). With these assumptions and caveats in mind, we
here discuss further the astrophysical significance of some of the features
that we have noted in the velocity maps in the previous section.

The maps shown in Fig. 22 can be directly compared to the corresponding
velocity maps from Gaia DR2 (Gaia Collaboration et al., 2018b, their figures
19 and 20). While many of the features out to about $3$ kpc were already
mapped with Gaia DR2, thanks to the larger sample of stars with line-of-sight
velocities, we can now map a significantly larger portion of the disc. The new
data reveals a complex and rich velocity field in the Galactic disc all the
way to the Galactic centre and the in-plane motions of the stars in the inner
Galaxy spectacularly reveal the clear signature of the Galactic bar. In
particular, we confirm the existence of a quadrupole pattern in the radial
velocity field, as previously found by Bovy et al. (2019) and Queiroz et al.
(2021). Additionally, we now clearly see the imprints of the Galactic bar in
the azimuthal velocity field, and in the stellar velocity dispersions. With
the additional guidance provided from of a numerical simulation, we are able
to estimate the aspect angle of the bar and its corotation radius through a
simple bisymmetric model of the kinematic signature of the bar in the
azimuthal velocity field. We find that the bar has an angle of $19.5\pm
2.5\degr$ with respect to the Sun-GC direction, as measured within $1.5<R\leq
2.8$ kpc, and a corotation radius of $5.4\pm 0.5$ kpc. Then, using the
observed angular frequency, we deduce a lower limit of the bar pattern
frequency of $\Omega_{\rm bar}=38.1^{+2.6}_{-2}$ $\rm\,km\,s^{-1}$ kpc-1, and
an upper limit for the location of the outer Lindblad resonance of $9.7\pm
0.5$ kpc within the framework of the epicycle approximation. Bland-Hawthorn &
Gerhard (2016) compiled previous estimates of the corotation radius and gave
$\Omega_{\rm bar}=43\pm 9$ $\rm\,km\,s^{-1}$ kpc-1, and a corotation radius
range of $4.5-7$ kpc. Applying the Tremaine-Weinberg method to their kinematic
sample, Bovy et al. (2019) estimated a bar pattern speed of $41\pm 3$
$\rm\,km\,s^{-1}$ kpc-1 and a corotation radius of $5.5\pm 0.4$ kpc. Other
works found $37.5$ $\rm\,km\,s^{-1}$ kpc-1 (Clarke et al., 2019), $39\pm 3.5$
$\rm\,km\,s^{-1}$ kpc-1 (Portail et al., 2017a) and $37.5-40$
$\rm\,km\,s^{-1}$ kpc-1 (Li et al., 2022), which is based on hydrodynamical
simulations and made-to-measure models. Our findings are thus in good
agreement with other estimates of the bar’s fundamental parameters. As for the
bar orientation angle, Bland-Hawthorn & Gerhard (2016) compiled a range of
$28\degr-33\degr$, while other recent estimates range from $20\degr$ to
$28\degr$ (e.g. Wegg et al., 2015; Portail et al., 2017b; Bovy et al., 2019;
Queiroz et al., 2021). Our value is thus at the low end of the proposed range.

The outer disk has a large variety of features and streaming motions, which
likely contain fundamental clues on the dynamical nature of the non-
axisymmetric structures in the Milky Way. Monari et al. (2016) modelled the
impact of the bar and a two-armed quasi-static spiral pattern in the Galactic
disc, both simultaneously and separately. The comparison of the observed
kinematic maps of the RGB stars confirms that the influence of the Galactic
bar likely dominates even the outer portions of the disc out to at least
$R=11$kpc. Faure et al. (2014) carried out test particle simulations to
predict the global stellar response to spiral perturbations in the Galactic
disc, in the absence of an external excitation (such as due to an accreting
satellite). They integrate stellar orbits in a 2-arm Lin-Shu type spiral
potential (without a bar) and produce maps of mean galactocentric radial
velocity ($V_{R}$). They show (their figure 6) that inside the spiral arms’
corotation, in the region traced by the arm, the mean $V_{R}$ is negative (of
the order of -7 $\rm\,km\,s^{-1}$), i.e., stars have a bulk motion towards the
Galactic centre. Meanwhile, in the region between the arms, the stellar radial
motion is positive, i.e., stars have a bulk motion towards the anticentre.
Outside corotation, the pattern is reversed. Trends in $V_{R}$ induced by the
spiral arms are also shown in Antoja et al. (2016) and other works. In our
sample of RGB stars, a spiral feature might be visible in the $V_{R}$ maps,
apparent as an elongated systematic positive feature, located internally to
the dashed line in Fig. 22 (left panel), marking the minimum interarm density.

Other features are also seen in the $\Delta V$ maps that are worthy of note.
Some of them might be explained as extinction artefacts (i.e. those that are
elongated along the line-of-sight), while others are likely real. In
particular, we see an asymmetry in all three components in the outer disk
beyond $R=10-11$ kpc in the RGB sample that does not as yet have a clear
explanation. This feature is likely related to an asymmetry already noted in
Gaia Collaboration et al. (2018b), where the authors follow the asymmetry in
the azimuthal and vertical components, ($V_{\phi}$ ,$V_{z}$) plane, showing a
clear bimodality of stars. Stars concentrate mainly in two clumps, one with
negative $V_{z}$ at lower $V_{\phi}$ which is more prominent in the north, and
one with positive $V_{z}$ at higher $V_{\phi}$ more visible in the south. The
different proportions of the clumps of the bimodality at different Z seems to
be the cause of the asymmetry.

For the OB sample, one of the main results is that the velocity field of the
OB stars shows streaming motions that have a characteristic length similar to
the spiral arm density. This has never been shown with such detail in
2-dimensional maps for OB stars. One should nevertheless bear in mind that our
sample of OB stars is expected to trace the local (out to $\approx 2-2.5$ kpc
in heliocentric distance) motions of the gas, rather than the large-scale
features of the Galaxy. Indeed, OB stars are typically young, and therefore
are expected to inherit the motion of the gas from which they were recently
born. However, due to the relatively small region sampled by our OB sample,
the streaming motions in the OB stars really map only the Local Arm.

A comparison between the OB and RGB stars’ velocity fields reveals numerous
differences. The observed differences are not unexpected, since these two
samples trace, respectively, dynamically cold and hot stellar populations of
the Milky Way. In the OB stars we clearly see streaming motions that are
associated with the spiral structure of this population. In contrast, the
signature of the spiral arms in the RGB sample is not clearly evident and, if
present, is seen in the radial motions. In any case, if there is a signature
of streaming motions related to the spiral arms in the RGB sample it is
consistent with a two-armed structure, possibly driven by the bar.

###  6.2 Caveats and shortcomings

No method or procedure is perfect, and here we have adopted one that is not
without imperfections. In this section we confess all of our sins.

In section 5 we construct maps of the mean velocity field in galactocentric
cylindrical coordinates based on the derived velocities in the same
coordinates of the individual sources. These in turn were derived, in part,
from the individual measures of the proper motions and estimated distances.
While our methodology is relatively straight-forward it has several
shortcomings which should be addressed in the future. Beyond a distance of
about three kiloparsecs, the largest source of uncertainty in our velocity
maps comes from the distance uncertainties, yet our treatment of these is not
completely satisfactory.

Our distances are taken from CBJ2021, which uses a Bayesian approach to
estimate the distances from both astrometry and photometry, incorporating a
prior that includes both current astrophysical knowledge of stellar structure,
as well as some informed assumptions about Galactic structure. The probability
distribution of the individual distances is in general asymmetric, but we have
rendered them symmetric in order to apply a traditional approach to
propagating the uncertainties to galactocentric coordinates that implicitly
assume that the uncertainties are symmetric. As discussed in sec 3.3, this is
in part done out of necessity as we are not able to sample the probability
distribution function distance for each star.

We have assumed that there is no correlation between the distance and proper
motion uncertainties, which cannot be true in as much as these distances are
informed by the parallax. On the other hand, as pointed out in section 3,
these correlations are likely to be unimportant at larger distances where the
distances are more constrained by the photometry than the astrometry.

When propagating our uncertainties we did take into account the correlation
between the proper motion components, which would result in a correlation
between the two velocity components in galactocentric coordinates. However, we
did not take these correlations into account when estimating the mean velocity
field in section 5.

Most grievously, when estimating the mean velocity and velocity dispersion for
a given volume element, we only consider stars found in that small volume
element, according to their estimated distances. That is, though we’ve taken
into account the uncertainties of the distances when deriving the
uncertainties in the individual velocities, we’ve implicitly assumed that that
these distances are perfect when binning the stars into cells for estimating
the mean velocity field and velocity dispersions.

Finally, it should be remembered that the CBJ2021 distance estimates are based
on a prior that includes a bar with an assumed geometry. As this prior becomes
increasingly important with increasing distance, there is reason for concern
that our velocity field may be influenced by this prior. Apart from the
CBJ2021 photogeometric distances, we have attempted using different distance
estimates: CBJ2021 geometric distances show very similar results, Starhorse3
(Anders et al., 2022) distances also show a similar bisymmetric behaviour,
while GSP-Phot distances have a similar trend but at a wrong distance. As
mentioned above, all these distance estimators have a prior that includes a
bar. In any case, we stress that this potentially introduces a bias in the
derived velocity field that is not considered in our qualitative comparison
with simulated data.

##  7 Conclusions

Using Gaia DR3 we have mapped the kinematics of the stars over an extensive
area of the Milky Way’s disc. This has been made possible thanks to the new
line-of-sight (radial) velocities, as well as new and reliable astrophysical
parameters that have allowed us to differentiate between young (OB) and old
(RGB) stars. While our sample of OB stars with 6D space and velocity data is
much more limited, their kinematics are seen to be distinctly different than
the RGB stars, likely reflecting the complex motions of the gas from which
they were recently born. Our RGB sample has allowed us to map the kinematics
of the Galaxy over nearly a quarter of the disc, providing us a first clear
picture of the large-scale kinematic signature of the bar. In order to
interpret the features seen in the kinematic maps we have relied on simple
comparisons with a simulation of a barred galaxy.

In contrast to the bar signature, clear evidence of streaming motions
associated with spiral arms is much less evident and, if present, is
consistent with the two-armed structure that we see in the near-infrared. With
regards to the non-axisymmetric structure in 3D configuration space, we have
made no attempt to reconstruct the stellar density. Nevertheless, we are able
to map the young OB sample and young clusters out to 4 to 5 kiloparsecs from
the Sun, confirming previous findings that the Local arm has a length of at
least 8kpc. While weaker than the inner Sag-Car arm or the outer Perseus arm,
it seems its nomenclature is perhaps not so appropriate, and that we should
return to its original designation as the Orion arm (van de Hulst et al.,
1954).

This study should only be considered as a preliminary exploration of what Gaia
DR3 has to offer with regards to galactic structure. Our kinematic maps are
derived from velocities of individual stars based on individual distance
estimates. A more ideal approach would be a derivation of the velocity field
recognising that the velocities, like the distances, should be inferred
quantities rather than derived quantities as we have treated them here.
Indeed, with regards to galactic dynamics, the individual velocities are not
of direct interest to us at all. Instead, an appropriate model of the mean
velocity field should be adopted and incorporated in the prior, then the
relevant parameters of the model (like the bar orientation) adjusted to arrive
at the most likely set of parameter values given the data. Alternatively, one
could generate a suite of simulations of the kinematics of the Milky Way,
assuming different parameters, model the uncertainties to transform these to a
suite of mock catalogues, then find the one that best matches the Gaia data.
In either case, any ”fitting” should ideally be done against the measurements,
in data-space rather than in model-space, taking into account the selection
effects on the data.

Though we have only taken a first look at the treasure of data that Gaia DR3
has to offer, it is clear that the Gaia mission continues to fulfil its
promise to provide the information needed to eventually arrive at a complete
understanding of the dynamical state and processes at work in shaping our
Galaxy. Certainly we can expect further discoveries and insights from the
community in the future as it digests this latest census of the stars in our
Milky Way.

###### Acknowledgements.

RD would like to thank Robert Benjamin for helpful discussions on the
treacherous history of mapping the spiral structure of our Galaxy. This work
presents results from the European Space Agency (ESA) space mission Gaia. Gaia
data are being processed by the Gaia Data Processing and Analysis Consortium
(DPAC). Funding for the DPAC is provided by national institutions, in
particular the institutions participating in the Gaia MultiLateral Agreement
(MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia
archive website is https://archives.esac.esa.int/gaia. Full acknowledgements
are given in Appendix D. This work has used the following software products:
TOPCAT, STIL, and STILTS (Taylor, 2005, 2006); Matplotlib (Hunter, 2007);
IPython (Pérez & Granger, 2007); Astropy, a community-developed core Python
package for Astronomy (Astropy Collaboration et al., 2018).

## References

  * Abdurro’uf et al. (2022) Abdurro’uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35
  * Abolfathi et al. (2018) Abolfathi, B., Aguado, D. S., Aguilar, G., et al. 2018, ApJS, 235, 42
  * Allen & Santillan (1991) Allen, C. & Santillan, A. 1991, Rev. Mexicana Astron. Astrofis., 22, 255
  * Alves (2000) Alves, D. R. 2000, ApJ, 539, 732
  * Anders et al. (2022) Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91
  * Anderson et al. (2019) Anderson, L. D., Wenger, T. V., Armentrout, W. P., Balser, D. S., & Bania, T. M. 2019, ApJ, 871, 145
  * Anderson et al. (2013) Anderson, R. I., Eyer, L., & Mowlavi, N. 2013, MNRAS, 434, 2238
  * Anderson et al. (2016) Anderson, R. I., Saio, H., Ekström, S., Georgy, C., & Meynet, G. 2016, A&A, 591, A8
  * Andrae et al. (2022) Andrae, R., Fouesneau, M., Sordo, R., et al. 2022, A&A, accepted
  * Antoja et al. (2018) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360
  * Antoja et al. (2020) Antoja, T., Ramos, P., Mateu, C., et al. 2020, A&A, 635, L3
  * Antoja et al. (2016) Antoja, T., Roca-Fàbrega, S., de Bruijne, J., & Prusti, T. 2016, A&A, 589, A13
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Babusiaux et al. (2022) Babusiaux, C., Fabricius, C., Khanna, S., & Muraveva, T. 2022, A&A, accepted
  * Bailer-Jones (2015) Bailer-Jones, C. A. L. 2015, PASP, 127, 994
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611
  * Bennett & Bovy (2019) Bennett, M. & Bovy, J. 2019, MNRAS, 482, 1417
  * Blaauw et al. (1960) Blaauw, A., Gum, C. S., Pawsey, J. L., & Westerhout, G. 1960, MNRAS, 121, 123
  * Bland-Hawthorn & Gerhard (2016) Bland-Hawthorn, J. & Gerhard, O. 2016, ARA&A, 54, 529
  * Blomme et al. (2022) Blomme, R., Frémat, Y., Sartoretti, P., Guerrier, A., & et al. 2022, A&A, accepted
  * Bono et al. (2005) Bono, G., Marconi, M., Cassisi, S., et al. 2005, ApJ, 621, 966
  * Bovy et al. (2015) Bovy, J., Bird, J. C., García Pérez, A. E., et al. 2015, ApJ, 800, 83
  * Bovy et al. (2019) Bovy, J., Leung, H. W., Hunt, J. A. S., et al. 2019, MNRAS, 490, 4740
  * Brown (2021) Brown, A. G. A. 2021, ARA&A, 59 [arXiv:2102.11712]
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Cantat-Gaudin et al. (2018) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018, A&A, 618, A93
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Carlin et al. (2013) Carlin, J. L., DeLaunay, J., Newberg, H. J., et al. 2013, ApJ, 777, L5
  * Castro-Ginard et al. (2019) Castro-Ginard, A., Jordi, C., Luri, X., Cantat-Gaudin, T., & Balaguer-Núñez, L. 2019, A&A, 627, A35
  * Castro-Ginard et al. (2022) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2022, A&A, 661, A118
  * Churchwell et al. (2009) Churchwell, E., Babler, B. L., Meade, M. R., et al. 2009, PASP, 121, 213
  * Clarke et al. (2019) Clarke, J. P., Wegg, C., Gerhard, O., et al. 2019, MNRAS, 489, 3519
  * Clementini et al. (2016) Clementini, G., Ripepi, V., Leccia, S., et al. 2016, A&A, 595, A133
  * Clementini et al. (2019) Clementini, G., Ripepi, V., Molinaro, R., et al. 2019, A&A, 622, A60
  * Clementini et al. (2022) Clementini, G., V., V., Garofalo, A., et al. 2022, A&A, submitted
  * Creevey et al. (2022) Creevey, O., Sordo, R., Pailler, F., et al. 2022, A&A, accepted
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197
  * De Silva et al. (2015) De Silva, G. M., Freeman, K. C., Bland-Hawthorn, J., et al. 2015, MNRAS, 449, 2604
  * De Somma et al. (2021) De Somma, G., Marconi, M., Cassisi, S., et al. 2021, MNRAS, 508, 1473
  * Deepak & Reddy (2018) Deepak & Reddy, B. E. 2018, AJ, 156, 170
  * Do et al. (2019) Do, T., Hees, A., Ghez, A., et al. 2019, Science, 365, 664
  * Drimmel (2000) Drimmel, R. 2000, A&A, 358, L13
  * Drimmel et al. (2003) Drimmel, R., Cabrera-Lavers, A., & López-Corredoira, M. 2003, A&A, 409, 205
  * Drimmel & Poggio (2018) Drimmel, R. & Poggio, E. 2018, Research Notes of the AAS, 2, 210
  * Faure et al. (2014) Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564
  * Feast & Catchpole (1997) Feast, M. W. & Catchpole, R. M. 1997, MNRAS, 286, L1
  * Ferrers (1877) Ferrers, N. 1877, Q.J. Pure Appl. Math, 14, 1
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Fouesneau et al. (2022) Fouesneau, M., Frémat, Y., Andrae, R., Korn, A., & et al. 2022, A&A, submitted
  * Fragkoudi et al. (2019) Fragkoudi, F., Katz, D., Trick, W., et al. 2019, MNRAS, 488, 3324
  * Gaia Collaboration et al. (2021a) Gaia Collaboration, Antoja, T., McMillan, P. J., et al. 2021a, A&A, 649, A8
  * Gaia Collaboration et al. (2018a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018a, A&A, 616, A1
  * Gaia Collaboration et al. (2021b) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021b, A&A, 649, A1
  * Gaia Collaboration et al. (2022a) Gaia Collaboration, Creevey, O., Sarro, L., Lobel, A., & et al. 2022a, A&A, accepted
  * Gaia Collaboration et al. (2018b) Gaia Collaboration, Katz, D., Antoja, T., et al. 2018b, A&A, 616, A11
  * Gaia Collaboration et al. (2022b) Gaia Collaboration, Klioner, S. A., Lindegren, L., et al. 2022b, A&A, accepted, arXiv:2204.12574
  * Gaia Collaboration et al. (2021c) Gaia Collaboration, Luri, X., Chemin, L., et al. 2021c, A&A, 649, A7
  * Gaia Collaboration et al. (2022c) Gaia Collaboration, Recio-Blanco, A., Kordopatis, G., et al. 2022c, A&A, accepted
  * Gaia Collaboration et al. (2021d) Gaia Collaboration, Smart, R. L., Sarro, L. M., et al. 2021d, A&A, 649, A6
  * Gravity Collaboration et al. (2022) Gravity Collaboration, Abuter, R., Aimar, N., et al. 2022, A&A, 657, L12
  * Gualandris & Merritt (2008) Gualandris, A. & Merritt, D. 2008, ApJ, 678, 780
  * Gum et al. (1960) Gum, C. S., Kerr, F. J., & Westerhout, G. 1960, MNRAS, 121, 132
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85
  * Hou (2021) Hou, L. G. 2021, Frontiers in Astronomy and Space Sciences, 8, 103
  * Hunt & Bovy (2018) Hunt, J. A. S. & Bovy, J. 2018, MNRAS, 477, 3945
  * Hunt et al. (2019) Hunt, J. A. S., Bub, M. W., Bovy, J., et al. 2019, MNRAS, 490, 1026
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90
  * Inno et al. (2021) Inno, L., Rix, H.-W., Stanek, K. Z., et al. 2021, ApJ, 914, 127
  * Jönsson et al. (2020) Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120
  * Jurić et al. (2008) Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864
  * Katz et al. (2019) Katz, D., Sartoretti, P., Cropper, M., et al. 2019, A&A, 622, A205
  * Katz et al. (2022) Katz, D., Sartoretti, P., Guerrier, A., & Panuzzo, P. 2022, A&A, submitted
  * Khanna et al. (2019a) Khanna, S., Sharma, S., Bland-Hawthorn, J., et al. 2019a, Monthly Notices of the Royal Astronomical Society, 482, 4215
  * Khanna et al. (2019) Khanna, S., Sharma, S., Tepper-Garcia, T., et al. 2019, MNRAS, 489, 4962
  * Khoperskov & Gerhard (2021) Khoperskov, S. & Gerhard, O. 2021, arXiv e-prints, arXiv:2111.15211
  * Laporte et al. (2020) Laporte, C. F. P., Famaey, B., Monari, G., et al. 2020, A&A, 643, L3
  * Leung et al. (2022) Leung, H. W., Bovy, J., Mackereth, J. T., et al. 2022, arXiv e-prints, arXiv:2204.12551
  * Levine et al. (2006) Levine, E. S., Blitz, L., & Heiles, C. 2006, Science, 312, 1773
  * Li et al. (2022) Li, Z., Shen, J., Gerhard, O., & Clarke, J. P. 2022, ApJ, 925, 71
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * López-Corredoira & Sylos Labini (2019) López-Corredoira, M. & Sylos Labini, F. 2019, A&A, 621, A48
  * Luri et al. (2018) Luri, X., Brown, A. G. A., Sarro, L. M., et al. 2018, A&A, 616, A9
  * Madore (1982) Madore, B. F. 1982, ApJ, 253, 575
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94
  * Marasco et al. (2017) Marasco, A., Fraternali, F., van der Hulst, J. M., & Oosterloo, T. 2017, A&A, 607, A106
  * Mathur et al. (2017) Mathur, S., Huber, D., Batalha, N. M., et al. 2017, ApJS, 229, 30
  * Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  * Monari et al. (2016) Monari, G., Famaey, B., Siebert, A., et al. 2016, MNRAS, 461, 3835
  * Monari et al. (2019) Monari, G., Famaey, B., Siebert, A., Wegg, C., & Gerhard, O. 2019, A&A, 626, A41
  * Morgan et al. (1953) Morgan, W. W., Whitford, A. E., & Code, A. D. 1953, ApJ, 118, 318
  * Peek et al. (2022) Peek, J. E. G., Tchernyshyov, K., & Miville-Deschenes, M.-A. 2022, ApJ, 925, 201
  * Pérez & Granger (2007) Pérez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21
  * Pietrukowicz et al. (2021) Pietrukowicz, P., Soszyński, I., & Udalski, A. 2021, Acta Astron., 71, 205
  * Poggio et al. (2020) Poggio, E., Drimmel, R., Andrae, R., et al. 2020, Nature Astronomy, 4, 590
  * Poggio et al. (2021) Poggio, E., Drimmel, R., Cantat-Gaudin, T., et al. 2021, A&A, 651, A104
  * Poggio et al. (2018) Poggio, E., Drimmel, R., Lattanzi, M. G., et al. 2018, MNRAS, 481, L21
  * Portail et al. (2017a) Portail, M., Gerhard, O., Wegg, C., & Ness, M. 2017a, MNRAS, 465, 1621
  * Portail et al. (2017b) Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017b, MNRAS, 470, 1233
  * Queiroz et al. (2021) Queiroz, A. B. A., Chiappini, C., Perez-Villegas, A., et al. 2021, A&A, 656, A156
  * Ramos et al. (2018) Ramos, P., Antoja, T., & Figueras, F. 2018, A&A, 619, A72
  * Reid & Brunthaler (2020) Reid, M. J. & Brunthaler, A. 2020, ApJ, 892, 39
  * Reid et al. (2019) Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3
  * Ripepi et al. (2022a) Ripepi, V., Catanzaro, G., Clementini, G., et al. 2022a, A&A, 659, A167
  * Ripepi et al. (2022b) Ripepi, V., Clementini, G., Molinaro, S., & Leccia, S. 2022b, A&A, submitted
  * Ripepi et al. (2019) Ripepi, V., Molinaro, R., Musella, I., et al. 2019, A&A, 625, A14
  * Rix et al. (2021) Rix, H.-W., Hogg, D. W., Boubert, D., et al. 2021, AJ, 162, 142
  * Romero-Gómez et al. (2015) Romero-Gómez, M., Figueras, F., Antoja, T., Abedi, H., & Aguilar, L. 2015, MNRAS, 447, 218
  * Romero-Gómez et al. (2019) Romero-Gómez, M., Mateu, C., Aguilar, L., Figueras, F., & Castro-Ginard, A. 2019, A&A, 627, A150
  * Rybizki et al. (2022) Rybizki, J., Green, G. M., Rix, H.-W., et al. 2022, MNRAS, 510, 2597
  * Sartoretti et al. (2022) Sartoretti, P., Blomme, R., David, M., & Seabroke, G. 2022, Gaia DR3 documentation Chapter 6: Spectroscopy, Gaia DR3 documentation
  * Spekkens & Sellwood (2007) Spekkens, K. & Sellwood, J. A. 2007, ApJ, 664, 204
  * Steinmetz et al. (2020) Steinmetz, M., Matijevič, G., Enke, H., et al. 2020, AJ, 160, 82
  * Tarricq et al. (2022) Tarricq, Y., Soubiran, C., Casamiquela, L., et al. 2022, A&A, 659, A59
  * Taylor & Cordes (1993) Taylor, J. H. & Cordes, J. M. 1993, ApJ, 411, 674
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Taylor (2006) Taylor, M. B. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 351, Astronomical Data Analysis Software and Systems XV, ed. C. Gabriel, C. Arviset, D. Ponz, & S. Enrique, 666
  * Tepper-Garcia et al. (2021) Tepper-Garcia, T., Bland-Hawthorn, J., Vasiliev, E., et al. 2021, arXiv e-prints, arXiv:2111.05466
  * Trick et al. (2019) Trick, W. H., Coronado, J., & Rix, H.-W. 2019, MNRAS, 484, 3291
  * Trick et al. (2021) Trick, W. H., Fragkoudi, F., Hunt, J. A. S., Mackereth, J. T., & White, S. D. M. 2021, MNRAS, 500, 2645
  * van de Hulst et al. (1954) van de Hulst, H. C., Muller, C. A., & Oort, J. H. 1954, Bull. Astron. Inst. Netherlands, 12, 117
  * Vasiliev & Baumgardt (2021) Vasiliev, E. & Baumgardt, H. 2021, MNRAS, 505, 5978
  * Wegg et al. (2015) Wegg, C., Gerhard, O., & Portail, M. 2015, MNRAS, 450, 4050
  * Widmark (2019) Widmark, A. 2019, A&A, 623, A30
  * Widrow et al. (2014) Widrow, L. M., Barber, J., Chequers, M. H., & Cheng, E. 2014, MNRAS, 440, 1971
  * Widrow et al. (2012) Widrow, L. M., Gardner, S., Yanny, B., Dodelson, S., & Chen, H.-Y. 2012, ApJ, 750, L41
  * Williams et al. (2013) Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101
  * Xiang et al. (2021) Xiang, M., Rix, H.-W., Ting, Y.-S., et al. 2021, arXiv e-prints, arXiv:2108.02878
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, MNRAS, 468, 3289
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579
  * Zari et al. (2021) Zari, E., Rix, H. W., Frankel, N., et al. 2021, A&A, 650, A112
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723

##  Appendix A Queries used to select samples

For the OB sample, the following query was used:

⬇

1SELECT g.*, ap.*

2FROM gaiadr3.gaia_source AS g

3INNER JOIN gaiadr3.astrophysical_parameters AS ap

4ON g.source_id = ap.source_id

5WHERE ((ap.teff_gspphot > 10000 AND (ap.spectraltype_esphs = ’O’ OR
ap.spectraltype_esphs = ’B’ OR ap.spectraltype_esphs = ’A’) AND ap.teff_esphs
IS NULL) OR (ap.teff_esphs > 10000 AND ap.teff_gspphot > 8000) OR
(ap.teff_esphs > 10000 AND ap.teff_esphs < 50000 AND ap.teff_gspphot IS NULL))

6AND power(g.parallax/100.,5) < power(10.,2.-g.phot_g_mean_mag+1.8*g.bp_rp)

For the RGB sample, the following query was used:

⬇

1SELECT g.*

2FROM gaiadr3.gaia_source AS g

3WHERE (g.teff_gspphot<5500 and g.teff_gspphot>3000) and (g.logg_gspphot<3.)

Subsequent cross-matches to the tables containing the CBJ2021 distances and
the astrometric fidelity flag were done after saving the results of the above
query as a user table. The final selection on $|Z|<300$ pc and $|Z|<1$ kpc for
the OB sample and the RGB samples, respectively, was done during the creation
of the maps in Section 4 and 5.

As written, these queries will pull all columns from the tables gaia_source
and astrophysical_parameters. To save space the user is recommended to pull
only the columns of interest.

##  Appendix B Maps of velocity uncertainties

###  B.1 Velocity uncertainties with Fisher formalism

We briefly explain how we estimate the uncertainties of the velocity $V_{k}$
from the negative log-likelihood given in Eq. (19): First, we minimise the
negative log-likelihood in order to obtain estimates of the mean velocity
$V_{k}$ and the dispersion $\sigma_{k}^{*}$. Second, we employ the Fisher
formalism by computing the Hessian matrix of second derivatives of Eq. (19)
w.r.t. the two fit parameters $V_{k}$ and $\sigma_{k}^{*}$ and the covariance
matrix of the two fit parameters is then given by the inverse of the Hessian.

We compute the first derivatives of Eq. (19) w.r.t. $V_{k}$ and
${\sigma_{k}^{*}}$, respectively, equate to zero and try to solve:

$\frac{\partial\mathcal{L}}{\partial
V_{k}}=-\sum^{N_{*}}_{j}\frac{v_{k,j}-V_{k}}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}=0\quad\Leftrightarrow\quad
V_{k}=\frac{\sum^{N_{*}}_{j}\frac{v_{k,j}}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}}{\sum^{N_{*}}_{j}\frac{1}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}}$
(22)
$\frac{\partial\mathcal{L}}{\partial{\sigma_{k}^{*}}}=\frac{1}{2}\sum^{N_{*}}_{j}\left[\frac{2{\sigma_{k}^{*}}}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}-2{\sigma_{k}^{*}}\frac{(v_{k,j}-V_{k})^{2}}{({\sigma_{k}^{*}}^{2}+{\sigma_{k}^{*}}_{v_{k,j}}^{2})^{2}}\right]=0$
(23)

Since ${\sigma_{k}^{*}}=0$ is not an allowed solution because the dispersion
has to be strictly positive, we can cancel out from this and obtain:

$\sum^{N_{*}}_{j}\frac{1}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}=\sum^{N_{*}}_{j}\frac{(v_{k,j}-V_{k})^{2}}{({\sigma_{k}^{*}}^{2}+{\sigma_{k}^{*}}_{v_{k,j}}^{2})^{2}}$
(24)

This cannot be solved analytically for ${\sigma_{k}^{*}}$ and we employ a
numerical solution instead.

Next, we need to compute the Hessian matrix of second derivatives of
$\mathcal{L}$:

$\frac{\partial^{2}\mathcal{L}}{\partial
V_{k}^{2}}=\sum^{N_{*}}_{j}\frac{1}{{\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}}$
(25) $\frac{\partial^{2}\mathcal{L}}{\partial
V_{k}\partial{\sigma_{k}^{*}}}=2{\sigma_{k}^{*}}\sum^{N_{*}}_{j}\frac{v_{k,j}-V_{k}}{({\sigma_{k}^{*}}^{2}+{\sigma_{k}^{*}}_{v_{k,j}}^{2})^{2}}$
(26)
$\frac{\partial^{2}\mathcal{L}}{\partial{\sigma_{k}^{*}}^{2}}=\sum^{N_{*}}_{j}\left[4{\sigma_{k}^{*}}^{2}\frac{(v_{k,j}-V_{k})^{2}}{({\sigma_{k}^{*}}^{2}+{\sigma_{k}^{*}}_{v_{k,j}}^{2})^{3}}-\frac{2{\sigma_{k}^{*}}^{2}}{({\sigma_{k}^{*}}^{2}+{\sigma_{k}^{*}}_{v_{k,j}}^{2})^{2}}\right]$
(27)

The last equation has already been simplified by using Eq. (24). We now write
all the second derivatives into the Hessian matrix:

$H=\left(\begin{array}[]{cc}\frac{\partial^{2}\mathcal{L}}{\partial
V_{k}^{2}}&\frac{\partial^{2}\mathcal{L}}{\partial
V_{k}\partial{\sigma_{k}^{*}}}\\\ \frac{\partial^{2}\mathcal{L}}{\partial
V_{k}\partial{\sigma_{k}^{*}}}&\frac{\partial^{2}\mathcal{L}}{\partial{\sigma_{k}^{*}}^{2}}\end{array}\right)$
(28)

Given the numerical best-fit solutions for $V_{k}$ and ${\sigma_{k}^{*}}$, the
Hessian $H$ is just a 2$\times$2 matrix of numerical values that we can
invert. According to the Fisher information criterion, the covariance matrix
of the fit parameters is equal to the negative of the inverse Hessian. Since
Eq. (19) is the negative log-likelihood, the covariance matrix is the inverse
Hessian without the minus sign:

$C=H^{-1}$ (29)

The standard deviation on $V_{k}$ is then simply given by $\sqrt{C_{00}}$.
Note that we cannot give an analytic expression for it because we cannot solve
Eq. (24) analytically for ${\sigma_{k}^{*}}$. Yet, we can compute
$\sqrt{C_{00}}$ from the numerical solution for ${\sigma_{k}^{*}}$.

###  B.2 Solution for ${\sigma_{k}^{*}}\gg\sigma_{v_{k,j}}$

The dispersion ${\sigma_{k}^{*}}$ may be much larger than all individual
measurement errors $\sigma_{v_{k,j}}$. In that case, we can approximate
${\sigma_{k}^{*}}^{2}+\sigma_{v_{k,j}}^{2}\approx{\sigma_{k}^{*}}^{2}$ such
that Eq. (22) and Eq. (24) simplify to:

$V_{k}=\frac{1}{N_{*}}\sum^{N_{*}}_{j}v_{k,j}\quad\textrm{and}\quad{\sigma_{k}^{*}}^{2}=\frac{1}{N_{*}}\sum^{N_{*}}_{j}(v_{k,j}-V_{k})^{2}$
(30)

In this specific case, there now is an analytic solution for $\sigma_{k}^{*}$.
The Hessian from Eq. (28) reduces to:

$H=\left(\begin{array}[]{cc}\frac{N_{*}}{{\sigma_{k}^{*}}^{2}}&\frac{2}{{\sigma_{k}^{*}}^{3}}\sum^{N_{*}}_{j}(v_{k,j}-V_{k})\\\
\frac{2}{{\sigma_{k}^{*}}^{3}}\sum^{N_{*}}_{j}(v_{k,j}-V_{k})&4\frac{\sum^{N_{*}}_{j}(v_{k,j}-V_{k})^{2}}{{\sigma_{k}^{*}}^{4}}-\frac{2N_{*}}{{\sigma_{k}^{*}}^{2}}\end{array}\right)$
(31)

Given $V_{k}=\frac{1}{N_{*}}\sum^{N_{*}}_{j}v_{k,j}$ we have
$\sum^{N_{*}}_{j}(v_{k,j}-V_{k})=0$ such that the off-diagonal elements of $H$
disappear and we are left with:

$H=\left(\begin{array}[]{cc}\frac{N_{*}}{{\sigma_{k}^{*}}^{2}}&0\\\
0&4\frac{\sum^{N_{*}}_{j}(v_{k,j}-V_{k})^{2}}{{\sigma_{k}^{*}}^{4}}-\frac{2N_{*}}{{\sigma_{k}^{*}}^{2}}\end{array}\right)=\left(\begin{array}[]{cc}\frac{N_{*}}{{\sigma_{k}^{*}}^{2}}&0\\\
0&\frac{2N_{*}}{{\sigma_{k}^{*}}^{2}}\end{array}\right)$ (32)

Since this Hessian has format 2$\times$2, we can easily invert it by hand to
obtain the covariance matrix $C$. In fact, we are only interested in the first
element of the diagonal of $C$:

$C_{00}=\frac{H_{11}}{H_{00}H_{11}-H_{01}H_{10}}=\frac{\frac{2N_{*}}{{\sigma_{k}^{*}}^{2}}}{\frac{N_{*}}{{\sigma_{k}^{*}}^{2}}\frac{2N_{*}}{{\sigma_{k}^{*}}^{2}}-0}=\frac{{\sigma_{k}^{*}}^{2}}{N_{*}}$
(33)

Consequently, the uncertainty on the mean velocity $V_{k}$ is approximately
$\frac{{\sigma_{k}^{*}}}{\sqrt{N_{*}}}$ if the dispersion satisfies
${\sigma_{k}^{*}}\gg\sigma_{v_{k,j}}$.

Figure 24: Velocity uncertainty maps for Milky Way RGB giants (left column)
and OB stars (right column). From top to bottom, the radial, azimuthal, and
vertical uncertainties are shown. The position of the Sun is marked by a
filled circle, and that of the Galactic Centre by a plus sign. The velocity
ranges have been chosen to enhance contrasts. These uncertainty maps can be
found in CDS in FITS format.

##  Appendix C Test particle simulations and Gaia mock catalogues

We use test particle simulations to guide the reader interpret the velocity
maps shown by the RGB sample in the inner regions and no dynamical modelling
has been undertaken. The initial conditions, the Galactic potential and the
steps performed in the integration process are described in Romero-Gómez et
al. (2015). Here we make a short summary.

Initial conditions for positions and velocities are drawn for a disc density
distribution following a Miyamoto-Nagai disc potential (Miyamoto & Nagai 1975)
with a typical scale-height ($h_{z}=300$ pc) and radial velocity dispersion
($\sigma_{U}=30.3\,$$\rm\,km\,s^{-1}$) of a red clump star. We first integrate
the initial conditions in the axisymmetric potential of Allen & Santillan
(1991) for $10\,$Gyr, we introduce the Galactic bar potential adiabatically
during $4$ bar rotations and we integrate another $4$ bar rotations so that
the particles get in statistical equilibrium with the final bar potential. The
Galactic bar consists of the superposition of two aligned Ferrers ellipsoids
(Ferrers 1877), one modelling the triaxial bulge with semi-major axis of
$3.13\,$kpc and the second modelling the long thin bar with semi-major axis of
$4.5\,$kpc. The bar rotates as a rigid body with a constant pattern speed of
$45\,$$\rm\,km\,s^{-1}$ kpc-1, placing corotation and the Outer Lindblad
Resonances at $4.7\,$kpc and $8.4\,$kpc, respectively.

In order to mimic the RGB sample, we generate two sets. One assuming for each
particle an absolute magnitude of $M_{K}=-1.61$ mag (Alves 2000) and no
dispersion, and another one with an absolute magnitude of $M_{K}=-3.0$ mag and
a dispersion of $1$ mag, corresponding to the Red Clump and the bright Red
Giant Branch, respectively. We assign an observed colour typical of a red
clump star (Alves 2000). We balance the particles in each set so that it
mimics the RGB proportion of Red Clump stars and brighter sources and the
apparent magnitude and spatial distribution. Using absorption (Drimmel et al.
2003) and extinction models (Cardelli et al. 1989), we obtain the Gaia G
magnitude. We follow the prescriptions described in the Gaia Science
performance webpage to obtain the astrometric and spectroscopic uncertainties.
In this case, since we use “photogeo” distances, we adopt a constant relative
uncertainty in distance of $15\%$, but we also made tests with a relative
uncertainty of $10\%$ and $20\%$.

The spatial distribution (top), and radial (middle) and azimuthal (bottom)
velocity maps corresponding to a density distribution with a Galactic bar
oriented $20^{\circ}$ from the Sun – Galactic centre line is shown in Fig. 25,
without (left) and with (right) uncertainties. As in the RGB sample, we only
show cells with a minimum of 20 stars, and mask all others. The solid black
line shows the orientation of the bar major axis in the simulation.

Figure 25: Spatial distribution (top), Median $V_{r}$ (middle) and median
$V_{\phi}$ (bottom) from a barred test particle simulation. Left: variables
without uncertainties. Right: variables with uncertainties. Dashed black
circle shows the corotation radius, black solid line the bar orientation angle
of $20^{\circ}$ and the black cross, the position of the Galactic Centre.

We fit the asymmetric velocity model of Eq. 20 to the mock velocity fields in
order to assess to which extent one can constrain the strength and phase angle
of the bisymmetric perturbation. This simple model is motivated by the
presence of periodic features seen in both $V_{R}$ and $V_{\phi}$ maps at low
radius, as caused by the dynamics of the bar. Furthermore, comparable forms
are regularly used to model the kinematics in barred galaxies (see e.g.
Spekkens & Sellwood 2007), although the axisymmetric radial component
$\overline{V}_{R}$ is rarely derived in the literature. For a disc in relative
equilibrium, this component should indeed cancel out as a function of radius.
However, with the current Gaia dataset is evident that this is not the case in
the azimuthal range covered by the RGB sample. This is the reason why we
decided to fit $\overline{V}_{R}$ as well with both mock and real kinematics.

As seen in Fig. 25, the angle $\phi_{\rm b}$ that the main axis of the
bisymmetry does with $Y=0$, the pseudo direction Sun $-$ Galactic centre, can
be measured independently using $V_{R}$ or $V_{\phi}$. Of course, this angle
is that of the bar at low radius. With $V_{R}$, this angle is the direction
where the right term of Eq. 20 changes its sign, that is $\phi_{\rm
b}=\phi_{R}-\pi/4$. With $V_{\phi}$ it is the angle where this term is
minimum, that is $\phi_{\rm b}=\phi_{\phi}-\pi/2$. An interesting feature seen
in the middle right panel of Fig 25 (mock data with uncertainties) is that the
direction of the change of sign of $V_{R}$ is no more perfectly aligned with
the bar major axis, but is closer to the Sun $-$ Galactic centre line. We thus
expect from the derivation of the simple asymmetric model that $\phi_{\rm
bisym}$ is smaller than the true bar orientation.

Figure 26 shows the amplitudes and phase angles of the bisymmetric model of
$V_{R}$ and $V_{\phi}$ applied to the simulation with a bar phase angle of
$-20\degr$, where we mask bins with $X>0$ to have the spatial distribution
almost similar to that of the Gaia data. The bottom and top panels
respectively show the results obtained with and without uncertainties added to
the mock variables. The asymmetric model applied to $V_{\phi}$ recovers well
the bar orientation within $R=2.2$ kpc, irrespective of the uncertainties,
while for $V_{R}$, the phase angle of the bar is only recovered when
uncertainties are not taken into account, and overestimates it in the other
case.