Gaia Data Release 3: Mapping the asymmetric disc of the Milky Way

With the most recent Gaia data release the number of sources with complete 6D phase space information (position and velocity) has increased to well over 33 million stars, while stellar astrophysical parameters are provided for more than 470 million sources, in addition to the identification of over 11 million variable stars. Using the astrophysical parameters and variability classifications provided in Gaia DR3, we select various stellar populations to explore and identify non-axisymmetric features in the disc of the Milky Way in both configuration and velocity space. Using more about 580 thousand sources identified as hot OB stars, together with 988 known open clusters younger than 100 million years, we map the spiral structure associated with star formation 4-5 kpc from the Sun. We select over 2800 Classical Cepheids younger than 200 million years, which show spiral features extending as far as 10 kpc from the Sun in the outer disc. We also identify more than 8.7 million sources on the red giant branch (RGB), of which 5.7 million have line-of-sight velocities, allowing the velocity field of the Milky Way to be mapped as far as 8 kpc from the Sun, including the inner disc. The spiral structure revealed by the young populations is consistent with recent results using Gaia EDR3 astrometry and source lists based on near infrared photometry, showing the Local (Orion) arm to be at least 8 kpc long, and an outer arm consistent with what is seen in HI surveys, which seems to be a continuation of the Perseus arm into the third quadrant. Meanwhile, the subset of RGB stars with velocities clearly reveals the large scale kinematic signature of the bar in the inner disc, as well as evidence of streaming motions in the outer disc that might be associated with spiral arms or bar resonances. (abridged)


Introduction
The determination of the structure and kinematics of the Milky Way has been investigated for more than a century. Researchers have been able to describe the morphology of external galaxies using deep photometric surveys, but the structure and evo-veys, such as the SDSS 1 , RAVE 2 , APOGEE 3 , LAMOST 4 , and GALAH 5 . Gaia has revolutionised this field, starting with first Gaia data release (DR1), providing new insights into the stability of the Galactic disc (e.g. Gaia Collaboration et al. 2018b; Antoja et al. 2018), its merger history (e.g. Helmi et al. 2018;Belokurov et al. 2018), and its structure through the discovery of new open clusters (e.g. Cantat-Gaudin et al. 2018;Castro-Ginard et al. 2019), to name a few. These results largely used the unprecedented number of about 7 million stars with 6D phasespace information in the second Gaia data release (DR2). The sample with Radial Velocity Spectrometer (RVS) measurements of Gaia Collaboration et al. (2018b), with full 6D phase-space measurements, already showed that the disc of the Milky Way is not kinematically axisymmetric. Since then, a large number of contributions have been published that provided new results and characteristics of the Milky Way disc. We refer to Brown (2021) for an updated review of the Milky Way with Gaia compared to a pre-Gaia view (e.g. Bland-Hawthorn & Gerhard 2016). The purpose of this contribution is to highlight the new information that is contained in the most recent Gaia data release regarding the structure of the Milky Way disc as revealed in configuration and velocity space.
Because it is rich in gas and has a disc structure, it was immediately expected that the Milky Way would have nonaxisymmetric structures, like other spiral galaxies. Clear evidence of spiral structure from the distribution of local OB associations was reported in the 1950s (Morgan et al. 1953). Since then, the location of spiral arms in density in the local neighbourhood has been studied using different tracers, such as giant molecular clouds, masers associated with high-mass star formation, H II regions, and young stars (OB stars, Cepheids, and young open clusters). Using parallaxes and proper motions of masers from the BESSEL survey, Reid et al. (2019) built logarithmic models of the Galactic spiral arms. The main arms identified in the model are the Norma-Outer arm, the Perseus arm, the Sagittarius-Carina arm (Sag-Car hereafter), and the Scutum-Centaurus arm. Another included arm is the Local Arm, which has been mostly considered as a minor feature with respect to the other arms listed above, as the name suggests. Maps showing the spiral arm segments of the Perseus, Local, and Sag-Car arms were produced by Poggio et al. (2021) (hereafter P21) with a local OB sample, Cepheids, and young open clusters, and similarly by Hou (2021) with spectroscopically confirmed OB stars. These maps extended the arm segments towards the third and fourth quadrant, where masers are mostly absent. While some progress has been made in detailing the large-scale spiral structure as evidenced by star formation products, the dynamical nature of these arms and the mechanisms causing their formation remains unknown.
The Galactic bar in the inner disc is another long-known asymmetry of the Galaxy, whose kinematic signatures can be found from the inner to the outer disc. Like many external barred galaxies, the Galactic bar has a boxy-bulge shape, but its length, orientation angle, and angular velocity are not yet well constrained. In this instance, we now have a strong asymmetry in the stellar distribution, at least for the inner regions of the Milky Way. While some evidence of asymmetry in the form of spi-ral arms extending farther out from the bar in the stellar disc is seen in the near-infrared (NIR) (Drimmel 2000;Churchwell et al. 2009), it has been challenging to find confirmation that the Milky Way hosts a density-wave-like structure in its kinematics. From earlier data releases, Gaia revealed that the velocity space of the stars is rich with structure. Most notably, the presence of arches and ridges in the V ϕ − R space (galactocentric azimuthal velocity and radius) indicate large-scale kinematic phenomena in the Galactic disc Ramos et al. 2018;Fragkoudi et al. 2019;Khanna et al. 2019). Disentangling these into identified resonances with the bar Fragkoudi et al. 2019;Trick et al. 2021;Monari et al. 2019;Laporte et al. 2020) and/or spiral arms and/or external perturbations is an ongoing process Khanna et al. 2019;Khoperskov & Gerhard 2021) and has proven to be difficult, mostly due to our detailed knowledge of the of stellar kinematics being contained within at most a few kiloparsecs of the Sun at most.
To assist us in understanding the kinematics of the Milky Way, comparisons between the observations and models will be important, and some of the previous works already mentioned have used this approach. Other recent works focused on how the spiral arms or a Galactic bar change the expected radial, tangential, and vertical kinematic maps (i.e. Faure et al. 2014;Monari et al. 2016;Hunt & Bovy 2018;Monari et al. 2019;Tepper-Garcia et al. 2021), based on either test particle simulations or pure N-body simulations. The observable used to compare with the data can either be directly mapping the average or dispersions of the velocity components on the Galactic plane, or checking the known moving groups in the solar neighbourhood or the ridges in the diagram of azimuthal velocity versus radius. In any case, an appropriate comparison of models to data must take the selection effects and uncertainties in the data into account. It is important to determine how they affect the prediction in contrast to ideal noise-free data. Romero-Gómez et al. (2015) showed the capabilities of the Gaia nominal mission in constraining the bar characteristics and constructing Gaia mock catalogues based on the Gaia science performance prescriptions for disc red clump stars.
In addition to the bar, spiral arms, and the Galactic warp, there has been some kinematic evidence of additional asymmetries that may indicate disequilibrium on a larger scale. With the RAVE survey, Williams et al. (2013) showed the presence of large-scale streaming motion in the disc and revealed differences above and below the Galactic plane. With SDSS data, Widrow et al. (2012) discovered similar wave-like compression or rarefaction features seen in both number density and bulk velocity, as well as towards the Galactic anticentre with LAMOST data (Carlin et al. 2013). The large-scale velocity field has also been mapped using highly precise line-of-sight velocities and distance tracers such as red clump giants (Bovy et al. 2015;Khanna et al. 2019a). (We adopt the term 'line-of-sight velocities' in place of 'radial velocities' for spectroscopically determined heliocentric velocity components in the direction of the source to avoid confusion with galactocentric 'radial velocity'.) Their results indicated streaming motion on scales much larger than about 2.5 kpc, but the analysis was likely limited by incompleteness in data coverage. As a demonstration of the enhanced astrometry and photometry in Gaia EDR3, Gaia Collaboration et al. (2021a) mapped the kinematics of the disc out to 14 kpc from the Galactic centre (GC). By selecting data in a narrow azimuthal range (20 • about the Galactic anticentre), they studied the azimuthal and vertical velocity components without requiring line-of-sight velocities. The large sample in their study allowed the dissection of the stellar rotation curve in the young and in the older pop-ulation of stars. Additionally, they showed that kinematic features (such as ridges in V ϕ − R space) seen in the inner disc with Gaia DR2 extended out to at least R = 14 kpc. By separately considering the stars above and below the Galactic plane, they also revealed that the lower disc has predominantly higher rotational velocities than the upper disc.
In this paper, we show the extraordinary capabilities of Gaia DR3 to shed light on the structure and kinematic issues mentioned above. We use similar tracers as in other works, using only the new information provided in Gaia DR3 to select our samples, and then to map the density and kinematics over a large portion of the disc. The paper is organised as follows. Section 2 describes the selection of the four tracers used in this work, namely clusters, Cepheids, OB stars, and red giant branch (RGB) stars, providing a description of their main properties. Section 3 describes the derivation of the positions, velocities, and uncertainties, including a short study of possible systematic effects. Section 4 maps the tracers into configuration space to show how they are distributed in relation to each other. Section 5 focuses on mapping the kinematics of the OB and RGB stars and on the information they contain about the bar and spiral arms. Section 6 discusses our results in context with other works and highlights the caveats and shortcomings that should be addressed in the future. In Section 7 we summarise our conclusions.

Selection of tracers
To map the asymmetry of the Galactic disc with Gaia, we selected young and old stellar populations: the former as the traditional tracer for the spiral structure that is used at optical wavelengths, where the surface brightness of disc galaxies, such as our Milky Way, are dominated by star formation products at optical wavelengths, while less-luminous older populations determine the mass distribution. The latest Gaia DR3 release allows us to select samples based on stellar parameters for the first time, which we used to select a sample of OB stars and red giants. The subset of sources with line-of-sight velocities  has full 6D phase-space information, allowing us to map the velocity field for these samples. In addition, we also investigate the distribution of open clusters and Classical Cepheids (DCEPs), for which we can derive excellent distances as well as ages. In this section we describe how we constructed each of these samples and how distances were derived for each.

Clusters
We used the list of probable members of the 2017 clusters studied by Cantat-Gaudin et al. (2020) with DR2 data. We obtained the DR3 source_id of these sources via the available crossmatch table, and removed stars whose EDR3 astrometry revealed them to be outliers by more than 3σ. This list of members was supplemented with the stars from 628 clusters that were recently discovered by Castro-Ginard et al. (2022) in the EDR3 catalogue, and the members found by Tarricq et al. (2022) in the outskirts of 389 nearby clusters. Most of these clusters have associated ages estimated with an artificial neural network applied to the Gaia DR2 data (for those in Cantat-Gaudin et al. 2020) or EDR3 photometry, and 988 of them are younger than 100 Myr.
The median astrometric parameters (parallax and proper motion) were computed for all clusters and are provided in Table 1. Before calculating the median parallax, we corrected the individual parallaxes following the recipe provided by Lindegren et al. (2021). Given the statistical precision obtained from using a large number of (corrected) parallaxes, we estimated distances 10 2 10 3 10 4 cluster distance [pc]  by inverting the median cluster parallax. The uncertainty on the bulk cluster astrometry was estimated as the quadratic sum of the statistical uncertainty and the uncertainty due to small-scale correlations, taken as 10 µas in parallax and 25 µas yr −1 in proper motion (Vasiliev & Baumgardt 2021).
We also computed the median line-of-sight velocity for the 2162 clusters in which at least one member had a DR3 line-ofsight velocity. Out of the 988 clusters younger than 100 Myr, 698 have line-of-sight velocities from DR3. The bulk line-ofsight velocities were computed from an average of 48 members per cluster, although this number varied significantly with age and distance (Fig. 1). For comparison, in Gaia DR2, the line-ofsight velocities were only available for an average of 10 stars per cluster. We estimated the line-of-sight velocity uncertainty as where σ v los,i are the nominal line-of-sight velocity uncertainties of the N cluster members, and 0.5 km s −1 is a conservative estimate of the line-of-sight velocity accuracy estimated in DR2 (Deepak & Reddy 2018; Katz et al. 2019). Although future investigations of the DR3 line-of-sight velocities are likely to show improved systematics with respect to DR2, this conservative choice has no significant impact on the results of this paper.

Classical Cepheids
The sample of DCEPs adopted in this work is mainly based on the list of sources in the vari_cepheid table, which is published in DR3 as a result of the processing by the Specific Objects Study (SOS) pipeline that was specifically designed to validate and fully characterise DCEPs and RR Lyrae stars observed by Gaia (hereafter referred to as SOS Cep&RRL pipeline) (see Clementini et al. 2016Clementini et al. , 2019Clementini et al. , 2022Ripepi et al. 2022b, for full details). This sample is composed of 3286 DCEPs belonging to the Milky Way, 1995 of which pulsate in the fundamental mode (F), 1097 in the first overtone (1O) and 194 are multimode (MULTI) pulsators. For these DCEPs, the SOS Cep&RRL pipeline provides pulsation periods, intensity-averaged magnitudes, peak-to-peak amplitudes, Fourier parameters, and other Notes. The full table is available as an electronic table via CDS. N: number of probable members kept to compute astrometric parameters. σ µα * , σ µ δ , and σ ϖ are the observed standard deviations of the members. σ v los : computed line-of-sight velocity uncertainty. N v los : number of members used to compute the cluster line-of-sight velocity. n 0 : number of members with an absolute magnitude brighter than 0.
quantities whose full description can be found in Ripepi et al. (2022b). The DR3 DCEPs sample was complemented with DCEPs taken by the recent compilations of Pietrukowicz et al. (2021) and Inno et al. (2021). We removed two additional multimode DCEPs pulsating in the second and third mode, as well as sources already in the DR3 DCEPs sample, and retained only objects with valid measurements of the mean magnitude in all three Gaia passbands and reliable proper motions. We find an additional 564 objects from Pietrukowicz et al. (2021) and 43 objects from Inno et al. (2021), with 27 objects in common.
For these 27, we adopted the classifications and periods from Pietrukowicz et al. (2021), giving a total of an additional 580 DCEPs from the literature. However, an additional 81 literature DCEPs were removed as suspect binaries from their position in the period-Wesenheit diagram (see next section). We were therefore left with 486 and 13 DCEPs from the Pietrukowicz et al. (2021) and Inno et al. (2021) catalogues, respectively. The total sample is therefore composed of 3785 DCEPs. However, as we describe below, we further clipped this sample.
2.2.1. Distances and cleaning of the sample.
An estimate of the distance to each DCEP in our sample was obtained directly from the definition of the distance modulus w− W=−5 + 5 log D, where w and W are the apparent and absolute Wesenheit magnitudes, respectively. The Wesenheit magnitudes are reddening free by construction, assuming that the extinction law is known (Madore 1982). The coefficient of the w magnitude has been derived in the Gaia bands on an empirical basis by Ripepi et al. (2019) and is defined as w = G −1.90×(G BP −G RP ). The absolute Wesenheit magnitude W was calculated using the period-Wesenheit-metallicity (PWZ) relation recently published by Ripepi et al. (2022a), To calculate the value of w for the DCEPs in our sample, we used different Gaia (G, G BP , G RP ) magnitude data sets. For the DCEPs in the DR3 vari_cepheid table, the SOS Cep&RRL pipeline provides intensity-averaged magnitudes in the three Gaia bands, that is, magnitudes that are calculated to best resemble the magnitude that the DCEPs would have if they were non-variable stars. Instead, for the 499 literature DCEPs that are not in the DR3 vari_cepheid table, we only have the mean magnitudes estimated in the Gaia photometric processing (see Riello et al. 2021, for details) that are available for all sources in the Gaia source catalogue. However, using mean magnitudes in the Gaia source catalogue for the literature DCEPs sample does not bias our results because it was found that the difference between w magnitudes calculated in the two different ways is only −0.01±0.03 mag (Ripepi et al. 2022a). Obtaining reliable values of w is possible only for sources with reliable values of the G, G BP , and G RP magnitudes. Objects with a magnitude close to or fainter than G = 20 mag are expected to have very poor G BP photometry, thus resulting in unreliable mean G BP magnitudes. We described how we cleaned the sample for this effect at the end of this section.
The other ingredient needed to calculate the distance to each DCEP is W, for which we need the period and iron abundance of each pulsator. The periods were taken from the vari_cepheid  Ripepi et al. 2022a, for details). Even if not particularly precise, the iron abundances obtained in this way allow us to use Eq. 2 to derive reliable distances. According to the PWZ relation in Eq. 2, the impact on the distance to a DCEP produced by an uncertainty of 0.11 dex in metallicity is ∼2.5%, and even considering a conservative uncertainty of 0.2 dex in [Fe/H], the uncertainty on the distance would be a still tolerable 5%.
After the values of w and W were estimated as explained above, it was straightforward to calculate the distance to each DCEP in our sample and its uncertainty. This was calculated by error propagation: σ d = 0.4605σ µ d, where σ µ is the uncertainty on the distance modulus, calculated by adding the uncertainties on w and W in quadrature, which in turn were estimated by propagating the uncertainties on the Gaia magnitudes for w, and for W from the uncertainties in the coefficients of Eq. 2, and the uncertainty on the iron abundance (the uncertainty on the periods is negligible).
An analysis of the derived distances and relative uncertainties revealed that many faint objects had unreliable distances and/or very large uncertainties as a result of the large photometric uncertainties, especially in the G BP band. To provide a cleaner sample for further analysis, we experimented with the data and reached the conclusion that retaining only DCEPs with distances smaller than 30 kpc and a relative distance uncertainty better than 10% is a good compromise between precision Fig. 2. Histogram of the G magnitude for the DCEP sample. 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. 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 that were useful for further analysis. The G-magnitude distribution of this selected sample is shown in Fig. 2. The number of faint objects, especially from the literature sample, has now been drastically reduced.

Line-of-sight velocities
The spectra collected with the RVS spectrometer on board Gaia allow measuring time-series of the line-of-sight velocity (RV) values for millions of stars with a G magnitude brighter than ∼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 DCEPs. (Fifteen and nine of these objects 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 verified 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±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. The uncertainties on these values can be evaluated on the basis of fig. 6.13 by Sartoretti et al. (2022).

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) devised a more accurate period-age-metallicity (PAZ) relation based on an updated theoretical pulsation scenario. Because the DCEP PWZ relation allows us to obtain individual accurate distances, it follows that the DCEP PAZ allows us to date any region in the Galactic disc in which a DCEP is present. To take advantage of this powerful tool, we adopted the following equation (see Table 9 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 wished to use the DCEPs to trace the Milky Way 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 properties of our selected DCEPs that are not published in other Gaia catalogues or papers.

OB stars
To select young stars on the upper main sequence, we used the effective temperatures provided in DR3, selecting stars with T eff > 10 000K. For hot stars, two sets of effective temperatures are provided in DR3 in general. One set is provided by a General Stellar Parameterizer from Photometry (hereafter GSP-Phot; see Andrae et al. 2022), which estimates stellar parameters using the Gaia G BP /G RP spectrophotometry, astrometry, and G band photometry. GSP-Phot makes different sets of parameter estimates using different stellar libraries, and for each source then chooses one of these as the best estimate. 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 uses the BP/RP spectrophotometry, without the astrometry, together with the RVS spectra if they are available as well (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. Only about half the sample of stars with T eff > 10 000K from either method has 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 used the spectral type determined by ESP-HS for all sources with Gaia BP/RP spectrophotometry as an additional assurance of quality. That is, T eff > 10 000K and the ESP-HS spectral type flag set to O, B, or A; For the stars with only ESP-HS temperatures, we required that the effective temperature be in the range 10 000 < T eff < 50 000K, as it was found that the small fraction of sources with T eff > 50 000K are likely to be unreliable ); For sources with both sets of stellar parameters, we required that the effective temperature T eff > 8000K for GSP-Phot and T eff > 10 000K for ESP-HS, letting the confirmation from GSP-Phot verify the sources with ESP-HS hotter than 50 000K. We note that we only used the effective temperature for the selection. A comparison of our temperatures against those found 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 σ logAge are the decimal logarithm of the age and its uncertainty.
in the literature (Mathur et al. 2017;Abolfathi et al. 2018;Xiang et al. 2021) for stars in the range 8000 to 10 000K shows an RMS difference smaller than 900K and offsets smaller than 400K. While the differences increase for higher measured T eff , 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 introduces undesired subdwarfs and white dwarfs in our sample, we imposed the additional criterion G +5 log(ϖ/100) < 2.+1.8(G BP −G RP ) to remove sources fainter than our target upper main-sequence stars. The colour term takes extinction and reddening into account, where 1.8 is approximately the slope of the reddening vector in (G, G BP − G RP ) space. To capture distant sources with negative parallaxes, we rewrote this criterion in the form (ϖ/100.) 5 < 10. (2.−G+1.8 * (G BP −G RP )) (4) (see Appendix A for an example query). These criteria together 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 removed these contaminants by keeping only stars whose distance from the Galactic plane is smaller than 300pc. This reduced the sample to 621 609 stars. Finally, we used 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 criterion removed only 7% of the sample, leaving us with a final sample of 579 577 stars, 91 836 (15.8%) of which have line-of-sight velocities. However, we note that the line-of-sight velocities of a fraction of them were estimated using an RVS spectral template with temperatures that are very different from the effective temperatures that were finally estimated for them, and therefore they are likely to have incorrect line-of-sight velocities . Removing stars whose rv_template_temp) is lower than 7000K gives us 77 659 stars with valid line-of-sight velocities. Figure 3 shows the G magnitude distribution of our OB sample.
For the purpose of mapping, we need distance estimates for our sources. While Gaia provides parallaxes, about 40% of our sample have significant (σ ϖ /ϖ > 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 ϖ/σ ϖ for our sample as well as for 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 astrometry and photometry. However, 43% of our sample are sources that have only temperature estimates from ESP-HS, and therefore they lack a distance estimate from GSP-Phot. We therefore adopt the photogeo distances from Bailer-Jones et al. (2021) that are based on astrometric and photometric data, as recommended in Fouesneau et al. (2022), and which are available for our entire sample 6 .
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 that are consistent with the mean OB star V ϕ curve derived below in section 5.

Giants
To select stars on the red giant branch (RGB), we used the effective temperatures and surface gravities provided in DR3, select-A&A proofs: manuscript no. output  ing stars with 3000 < T eff < 5 500K and logg < 3.0, as provided by GSP-Phot  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 or PHOENIX libraries ). 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.
As in the OB sample, and in order to perform density and velocity maps, we needed to choose a distance estimator. The distribution of ϖ/σ ϖ 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 ϖ/σ ϖ < 5. All stars in the giant sample have GSP-Phot parameters. One option can therefore be to use the provided distance in Gaia DR3 ). Another option are the geo and photogeo distances from Bailer- Jones et al. . In order to choose the appropriate distance estimator for the large extent of the RGB sample, we cross-matched 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. This resulted 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). As shown in Babusiaux et al. (2022), the large parallax uncertainties and the prior used moreover causes the derived GSP-Phot distances to concentrate in density, forming a ring around 2 kpc. This makes them inappropriate for studying the inner disc. Therefore, and as in the OB sample, we adopted the photogeo distance estimate for the RGB sample.
In addition, as in the selection of OB stars, we kept only sources with good astrometry, that is, with an astrometric fidelity f a > 0.5. This criterion removed only 14% of the sample, leaving us with a sample of 9 959 807 stars, 6 586 329 of which have lineof-sight velocities.
We found that the above selection also included many sources from the Large and Small Magellanic Clouds, as well as a number of globular clusters. Because the goal of the paper is to study the Galactic disc, we performed an additional cut on the altitude with respect to the Galactic plane, that is, |Z| < 1 kpc. We also explored whether this selection might 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 at the faint end of our sample and therefore are a negligible fraction in our full sample, and are not expected at all in the sample with RVS line-of-sight velocities, so we did not attempt to remove 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 we used throughout the paper. We provide the number of sources with full 5D astrometry and the subsample including 6D phase-space information (astrometry and line-ofsight velocities).

Mapping to configuration space
To map our tracers in a 3D Cartesian coordinate space, we must transform astrometric angular measurements (ϖ, α, δ) into 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 Fig. 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). 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 adopted the photogeo distances of CBJ2021, for the clusters, we took the inverse of their median parallax, and for the Cepheids, we used a photometric distance based on the Leavitt law, as detailed in Section 2.2.
It has been known for a long time that some Milky Way DCEPs are members of OCs (see e.g. Anderson et al. 2013, and references therein). It is therefore useful to compare the photometric distances inferred for DCEPs 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 crossmatched the DCEP 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σ in all the cases, except for the farthest OC of the sample, namely UBC 608. The discrepancy for this source is smaller than 1.5σ, however.
From the heliocentric distances d, we can easily derive heliocentric Cartesian coordinates (x, y, z) under the assumption of a (non-relativistic) 3D Euclidean geometry. Because we assumed the positions to be known, we took advantage of the provision of (l, b) of each tracer in the Gaia archive and used the usual transformations, where positive x is towards the Galactic centre. Galactocentric Cartesian coordinates are then typically derived as a simple translation of the origin, where R ⊙ and Z ⊙ are the distance of the Sun from the Galactic centre and above the Galactic (Z = 0) midplane.
For our choice of R ⊙ , we used 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 ⊙ = 8277 ± 9(stat) ±30(sys) pc, although we note that this value disagrees within the uncertainties with the independent determination of R ⊙ by Do et al. (2019) (7.959 +/-59 (stat) +/-32 (sys)) and also with their previous determinations, indicating that our assumed R ⊙ may in fact be in error by as much as 200 to 300 parsecs. Our adopted value of R ⊙ assumed 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 ⊙ = 8.23 ± 0.12kpc based on observed stellar kinematics towards the Galactic centre. This value 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 ⊙ only results in a trivial offset in the maps. However, as discussed further in section 3.4, it can introduce rather undesirable effects when the velocities are mapped in galactocentric cylindrical coordinates.
The transformation into galactocentric coordinates (Eq. 6) is an approximation because 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 Blaauw et al. 1960), there may well be a residual offset due to the height of the Sun above the Z = 0 plane. It was already noted at the time that determinations of Z ⊙ from hydrogen radio emission do not coincide with those based on nearby stellar  ) rather complicates this discussion because the local stellar mid-plane of the Galaxy might very well not coincide with a Z = 0 plane as defined by the average vertical density distribution of the inner disc. 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 ⊙ . As a case in point, using a very local sample of stars, Gaia Collaboration et al. (2021d) reported that Z ⊙ 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 δθ between the galactic (b = 0) plane and the Galactic (Z = 0) plane of about only 0.1 • (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 ⊙ = 0, so that Z = z.

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 contain 77 659 and 5 730 578 sources, 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).
The relative velocity components in the heliocentric Cartesian coordinates defined by Eq. 5 are where (µ α * , µ δ ) are the proper motion components, v r is the lineof-sight (i.e. radial) velocity, A ′ 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 star, Alternatively, the astrometry might first be converted into galactic coordinates by removing A ′ from Eq. 7 and substituting the proper motion components in ICRS with those in Galactic coordinates, and then computing v rel from the Galactic coordinates.
In this case, the triad (Eq. 8) would be in (l, b) rather than (α, δ). 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 los = radial_velocity − 7.98 + 1.135grvs_mag.
This was applied to stars for which 8500 ≤ rv_template_teff ≤ 14 500 K and 6 ≤ grvs_mag ≤ 12. The velocities v rel = (u, v, w) derived above are relative to the Sun. To place them in a galactocentric reference frame, we must add the solar velocity with respect to the Galactic centre, v ⊙ , Traditionally, the solar velocity v ⊙ 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, which commonly is assumed to be in circular motion about the Galactic centre. However, with the recent precise measurement of the proper motion of the SagA*, together with R ⊙ , the azimuthal and vertical components of the solar galactocentric velocity can be derived in a more direct and precise way. From Reid & Brunthaler (2020), we have (µ l , µ b ) = (−6.411 ± 0.008, −0.219 ± 0.007) mas yr −1 for the proper motion of SagA*, which together with R ⊙ gives the solar Y and Z-velocity components. The same reduction of the SagA* S2 data by the GRAVITY Collaboration that yielded R ⊙ also yields the line-of-sight velocity towards SagA*, interpreted as the reflex motion of the solar velocity towards the Galactic centre. This results in if we assume that SagA* is stationary with respect to the Galactic centre. (see , for further discussion on this approach to deriving v ⊙ ). The uncertainties in v ⊙ as well as any error in our adopted R ⊙ gives us a systematic error common to all our galactocentric velocities v * . See section 3.4 below for further discussion.
The (u, v, w) components of our galactocentric velocities v * are rigorously in the same coordinate system as defined by Eq. 5, that is, they are 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 δθ mentioned above, if Z ⊙ 0. However, the systematic error introduced by ignoring this (unknown) tilt is much smaller than the systematic error introduced by the uncertainties in v ⊙ .
Assuming Z ⊙ = 0, the galactocentric radial and azimuthal velocities can be found from where ϕ is the galactocentric azimuth, taken positive in the direction of Galactic rotation, making our (R, ϕ, z) galactocentric cylindrical coordinates a lefthanded system.

Propagation of uncertainties
To estimate the uncertainties in positions and velocities from the formal errors in astrometry, we chose to concatenate the Jacobian matrices of the consecutive transformations necessary to move from the initial reference frame to the desired one. This means that we implicitly linearised the functions that allow us to convert astrometry into positions and velocities. In other words, we simplified the sequence of non-linear transformations (e.g. see eqs. 7, 12, and 13) that convert the coordinates in one reference frame, x 1 , into the coordinates in another frame, x 2 , with a single matrix product of the form 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 into the next, where f i is the function that calculates the i th th component of the vector x 2 in the desired coordinate frame given the vector x 1 in the original coordinate frame.
In practice, we first constructed the covariance matrix in the frame of the Gaia astrometry using the proper motion and lineof-sight velocity uncertainties and their associated correlations. As said above, we neglected the uncertainties in sky position, and consequently, all involved terms were set to zero. At the same time, we replaced the uncertainty in parallax with the uncertainty on our distance estimate, σ d (discussed below), and set the corresponding correlations to zero. The resulting initial covariance matrix for any source is, thus, where ρ µ * α µ δ 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), We then again used equations 17 and 15 to successively move first into galactocentric Cartesian (eq. 6 and eq. 10) and, finally, into galactocentric cylindrical coordinates (eqs. 12 and 13). At each step, we obtained 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, 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 is also a source of correlations regardless of the coordinate frames. However, the correlations introduced by the coordinate transformations in concert with the distance uncertainties 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 mock observables based on the covariance matrix in the initial Gaia frame. However, this approach contains some shortcomings: i) it requires a large number 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 because 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 implicitly assumed 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 where r_hi_photogeo and r_lo_photogeo are the 84th and 16th percentiles, respectively, of the photogeo distances provided by CBJ2021. Figure 9 shows the variation in the distance uncertainties in our sample with distance for the subset of RGB stars with radial velocities. We note that the median relative uncertainties are smaller than 20%. The distance uncertainties for the OB sample with line-of-sight velocities follows the same trend, but covers 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. These correlations are not expected in purely photometric 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, which we used for the OB and RGB samples, the relatively nearby sources, which generally have small relative distance uncertainties, are constrained by the parallaxes, while the more distant sources with larger uncertainties are 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 considered only the correlations between the proper motion components in the initial covariance matrix (eq. 16).
We note that because the tangential (perpendicular to the line-of-sight) velocity components are dependent on both distance and proper motion, their uncertainties are perfectly correlated with the distance uncertainties. As the line-of-sight velocity is not correlated with the astrometry and may be of a different magnitude than the uncertainty in the tangential velocity, the correlations and uncertainties of our velocity components consequently 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, is dominated by the uncertainty in the line-of-sight velocities. On the other hand, the uncertainties for the two 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 into velocities in galactocentric cylindrical coordinate, via equations 12 and 13. Errors in ϕ can become quite large near the Galactic centre. Moreover, by performing the coordinate transformation to (v R , v ϕ ) above, we have made ourselves vulnerable to systematic error in our assumed values of R ⊙ and v ⊙ . We explore the effect of these types of uncertainties in the following section.

Effect of systematic errors
As mentioned above, while the formal errors are quite small, there may be significant systematic errors in our assumed values for R ⊙ and v ⊙ . In addition, as we demonstrate below, even random distance errors can introduce systematic errors in the mean galactocentric velocity components. To investigate the possible effects that these errors can introduce, we constructed a mock catalogue from a rather artificial distribution of stars: We uniformly populated a disc of stars centred on the Sun with a radius of 8 kpc, and added a Gaussian velocity dispersion and an azimuthal velocity with respect to the Galactic centre at R GC that was consistent with what we observe for the RGB sample (see figure 17), but assumed no mean radial motion (i.e. v R = 0). We also assumed a motion of the observer at the Sun and derived the proper motion and line-of-sight velocities from the relative motions, to which we added fractional uncertainties of 0.01 in the proper motions and 0.1 in the radial velocities. We then rederived the (v R , v ϕ ) velocity components from the observed (noisy) proper motions, distances, and line-of-sight velocities, and constructed maps of the observed mean velocity field, as described in section 5.
To model the effect of distance uncertainties, we added a 20% Gaussian uncertainty to the true parallax, and then took the inverse of the observed parallax as our distant estimate. We note that this is much larger than the actual uncertainties of our data set, which only reach a relative uncertainty of 20% of the distances at about 10 kpc from the Sun. It 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 towards larger distances. More importantly, for our rather artificial uniform disc of mock stars, the mean estimated distance will be slightly underestimated systematically 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 Observed V R assuming that the true value of R ⊙ is 200 pc higher than the adopted value, together with the distance uncertainty model explained in the text. Lower right panel: Observed V R assuming that the true value of the solar velocity component V ϕ,⊙ is 10 km/s higher than the adopted value, together with the distance uncertainty model explained in the text. inferred mean velocity field (V R , V ϕ ), as shown in the upper two panels of figure 11. In the following, we refer to these 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 σ ϖ /ϖ = 0.02d kpc (i.e. σ ϖ = 20µas), we investigated the effect of assuming erroneous values for R ⊙ and the velocity of the Sun with respect the values that were used to generate the mock observations. We first introduced a systematic error in our assumed value for R ⊙ , 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 transformed the observed proper motions, distances, and line-of-sight velocities into (v R , v ϕ ) 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 ⊙ consistent with the mock velocities, but assuming an incorrect velocity for the solar v ϕ component, a similar but inverted pattern in the inferred V R velocities is observed, 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 is 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 sharply decreases beyond some limiting distance. Beyond this distance the sources used to estimate velocities have systematically overestimated 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 as a result of 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 caution 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.

Coordinate maps
In this section, we map the spatial distribution of the OB and RGB stars, open clusters, and Cepheid 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 in which the stellar density is markedly higher than in others. These 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, but with higher angular frequencies, as the density here is not derived using a smoothing kernel. In addition to 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 because 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 towards the Galactic centre, presumably due to the Galactic bar. Additionally, we note that an overdense region is apparent at x ≃ 2-3 kpc, running across a range of y. This is due to the combination of two main factors: (i) we used a magnitude-limited sample, implying that the density decreases with heliocentric distance, and (ii) the intrinsic distribution of RGB stars in the Galactic disc increases towards the inner parts of the Galaxy (as expected from an exponential disc).
To better explore the Galactic spiral structure, we applied the same approach to our OB sample as was adopted by P21 to map the stellar overdensity, defined as where the local surface density Σ(X, Y) and the mean surface density ⟨ Σ(X, Y) ⟩ were 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. The red diagonal stripes in Figure 13 correspond to the segments of the nearest spiral arms, consistent with the features identified in the left panel of Figure 12. Based on Figures 12 and 13, the emerging picture of the local Galactic spiral structure agrees well with what was found in P21 (see their Fig. 1B and 1C). The same Figure 13 shows a comparison between the OB overdensity map and the young (< 63 Myr) and bright (i.e. with at least five 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. Figure 13 shows a good agreement between the open clusters distribution and the spiral structure mapped by the OB sample. Figure 14 shows the spatial distribution of clusters younger than log t=7.8 (∼63 Myr). Although the young clusters clearly trace multiple elongated structures, evocative of arms and interarm 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) very nicely traces the orientation of the Perseus arm in the upper main-sequence stars (see also P21), and appears to agree reasonably well with the distribution of the OCs as well (see the dashed orange line in Fig. 14). According to this model, the two groups of OCs that Reid et al. (2019) used 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.
Finally, we mapped 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, and P21). The right panel of Figure 15 shows a comparison between the Cepheids WT and some models available in the literature. Solid lines show the four-armed model of Taylor and Cordes (Taylor & Cordes 1993); in the region in which the Cepheids are more abundant (x < 5 kpc), the Sagittarius Carina arm appears to agree well below approximately y ≃2.5 kpc. At approximately x ≃ 2.5 kpc and y ≃ 5 kpc, the overdensity of Cepheids seems to diverge with respect to the Taylor and Cordes model of the Sag-Car arm in this direction. 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 with 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.

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 used the selections of OB and RGB stars within |z| ≤ 0.3 (77 659 stars) and |z| ≤ 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 summarises our notation and is provided for convenience.

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 × 341 pixels. 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 2D histogram using cells of 100x100 pc. The Galactic centre is towards the right and is shown with a cross in the right-hand plot. Galactic rotation is clockwise. The position of the Sun is shown with a black dot. A logarithmic stretch is used for the map of the RGB stars to enhance lower-density regions. 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 at least 20 stars, and masked all others. This resulted in grids with a median number and maximum number 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 to a cell j of 100 × 100 pc in our three velocity components k = R, ϕ, or z maps, which contains N * stars. We estimated the stellar mean velocity V k and its associated dispersion σ * k by optimising the

S c u t u m S a g i t t a r i u s L o c a l P e r s e u s C y g n u s
Fig. 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σ uncertainty on the distance. They take statistical and systematic parallax errors into account. log-likelihood of the distribution of observed velocities v k, j (with uncertainties σ vk, 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 Article number, page 15 of 36 A&A proofs: manuscript no. output mean azimuthal velocity at R σ * R , σ * ϕ , σ * z velocity dispersions in R, ϕ, z directions Appendix B.1 details the derivations of the uncertainty σ V k, j of our mean velocities V k, j , which also demonstrates that under certain conditions, we can approximate them as σ V k = σ * k / √ N * . (Most inconveniently, σ 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.) 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 visually identify regions where streaming motions occur. The left panels of Fig. B.1 show the associated uncertainty maps for the velocities.
The V R map shows a remarkable bisymmetric feature on either side of the GC, with negative and positive values on each side of the apparent major axis of the bar. This quadrupole feature is a characteristic of the mean inward motion down to ∼ −40 km s −1 (y > 0) and the mean outward motion up to ∼ 45 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 this quadrupole pattern in detail. We also note that the RGB sample contains enough stars to apparently provide us with mean velocity estimates beyond the Galactic centre (GC), but they should be interpreted with caution.
The σ * R map shows a bisymmetric pattern as well, but different from V R , with larger amplitudes that are aligned with the direction in which V R changes its sign along the major axis of the bar, 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, for example R = 6.5 kpc (x ∼ 1.5 kpc), which shows that V R is larger for |y| < 2 kpc than at smaller and larger azimuths. At R ∼ 10 kpc (x, y ∼ −1.5, +1 kpc), V R shows a clear change of sign with respect to azimuth.
The distribution of the azimuthal velocity V ϕ 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 σ * ϕ map also seems to exhibit a bisymmetry in the bar region, but rotated by about 45 • with respect to that seen in σ * 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 ≤ x ≤ 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 towards the Galactic anticentre, clearly showing the kinematic signature of the warp of the Galactic disc beyond, in agreement with previous results ( . We do not observe any bisymmetric feature in the vertical dispersion, which is larger along l ∼ 0. The vertical velocity component is more sensitive to systematic errors because our sources are near the galactic (b = 0) plane, therefore this component is predominantly seen 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 Milky Way 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.

Analysis of radial profiles
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 V ϕ (R) at that radius. We discarded radial bins with fewer than five pixels from the maps. 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 V ϕ ∼ 220 km s −1 to R ∼ 6 kpc and 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 their asymmetric drift is smaller, and it decreases with galactocentric radius. It rotates 17 km s −1 faster on average than the RGB stars over the radial range R = 6.5 − 10 kpc. These average rotation curves were subtracted from the V ϕ maps to produce residual velocity fields V ϕ (x, y) − V ϕ (R). These are useful to show 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 were calculated using overlapping radial bins of 1 kpc step and 2 kpc width. As expected, the agreement between the OB stars and the OCs is good. The Cepheids exhibit a mean azimuthal velocity similar to but slightly slower than 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 an asymmetric drift for an older population with a higher velocity dispersion.
The radial and azimuthal velocity dispersion profiles of RGB stars show two distinct regions. Within the inner part of the bar (R ∼ 2.5 kpc), the profiles are shallow and the dispersions comparable (∼ 75 − 80 km s −1 ), but we note that the (x, y) map of the dispersion (Fig. 16) in this region shows very strong asymmetries. Beyond 2.5 kpc, σ * z < σ * ϕ < σ * R and the profiles continuously decrease for all three components. The radial profiles of the velocity dispersion σ * 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 km s −1 , respectively. This agrees very well with the observed random motions of gas within R = 8 kpc (4−9 km s −1 ; Marasco et al. 2017). The axis ratios of the velocity ellipsoid as averaged from these profiles for R > 3 kpc are (σ * ϕ /σ * R , σ * z /σ * R , σ * z /σ * ϕ ) = (0.81, 0.60, 0.75) for the RGB stars and (0.66, 0.44, 0.66) for OB stars.

Kinematics of the bar
In this section, we estimate some fundamental parameters of the Galactic bar, guided by the 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 wish to clarify that this is not a made-to-measure, customised simulation that intends to reproduce the observed data set of RGB stars quantitatively. The simulation has to be considered as a simple diagnostic tool for a qualitative comparison to the observations, and the results presented below are tentative possibilities. We refer to Appendix C for the description and analysis of the simulation. We thus applied the same recipes as described in Appendix C to determine by analogy the orientation and pattern speed of the Galactic bar, as well as the location of the outer Lindblad resonance.
We performed a fit to 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 in the maps. This approximation of V R and V ϕ is referred to as V R,mod and V ϕ,mod , and is given by where V R (R) is the axisymmetric mean value, and A R and ϕ 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 ϕ,mod . 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 applies best to regions in which the kinematics is well described by a second-order perturbation, we restricted the analysis to 0 ≤ R ≤ 7 kpc and considered a radial bin width of 200 pc. Defining the residual velocity as V R,res = V R − V R,mod , the conditional likelihood function for each radial bin is expressed as and similarly for ϕ component, where ξ 2 = σ 2 V R + V 2 R,s , and n pix is the number of pixels inside the corresponding radial bin, requiring n pix > 25. Radial bins that did not satisfy this condition were 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 were uniform and spanned [−50, 50] km s −1 for V R , [0, 300] km s −1 for V ϕ , [0, 50] km s −1 for A R , [0, 30] km s −1 for V R,s , and [−80, 110] • for ϕ R . Following prescriptions from Appendix C, ϕ R and ϕ ϕ were linked to the direction of the bisymmetric perturbation of density with respect to the Sun-GC direction in the bar region, which we call ϕ b , by ϕ b,R = ϕ R − π/4 in the case of the V R model, and ϕ b,ϕ = ϕ ϕ − π/2 in the case of V ϕ . We quote the uncertainties on the parameters at the 16th and 84th percentiles of the posterior distributions. Figure 19 shows the resulting fits. 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 ∼ 0.8 and 6 kpc for A ϕ . As expected, this simple model shows that V R is anything but axisymmetric, as A R exceeds V R . The axisymmetric radial velocity is mostly positive, showing a bulk outward motion of 3.6 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 that is nearly relaxed, we attribute this non-negligible 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 parameters V R,s and V ϕ,s , which is ∼ 6 km s −1 at small radius and decreases to < 1 km s −1 out to R = 7 kpc (not shown in Fig. 19). As for V ϕ , we find that 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 ϕ models, particularly beyond R = 3 kpc, as shown in the bottom panel of Fig. 19. At the peak of the strength at R = 2.5 kpc, ϕ b differs by ∼ 12 • . Interestingly, ϕ b in the V ϕ model remarkably shows the same trend as that seen in the numerical simulation (Fig. C.2 of Appendix C), where it faithfully traces the true bar orientation even in the presence of uncertainties. As shown in Fig. C.3, the ϕ 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 ϕ b,ϕ remains flat before the abrupt drop, that is, −19.2 • ± 1.5 • for 1.5 < R ≤ 2.8 kpc. Then, as the location of the minimum of ϕ b of the V ϕ 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. C.2), we find by analogy that the range R = 5.2 − 6 kpc hosts the corotation radius of the Galactic bar, R CR . We adopted a conservative value R CR = 5.4 ± 0.5 kpc. This is surprisingly also the location at which ϕ b 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.
In Fig. 20 we show the angular velocity curve, Ω = V ϕ /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 was derived by assuming the epicycle approximation, which might not be perfectly correct in a radial range far from the solar neighbourhood. When we assume the corotation radius range obtained above from the phase of the bisymmetry in V ϕ , the intersection between R CR and the Ω curve provides the bar pattern speed, which is Ω bar = 38.1 +2.6 −2. km s −1 kpc −1 . We can also provide an estimation of the outer Lindblad resonance (OLR), which is the galactocentric radius at which the fixed pattern speed intersects the curve Ω + κ/2, that is, the position in the disc at which stars rotate slower than the bar pattern and perform two radial oscillations for one revolution of the bar pattern, κ being the epicyclic frequency. In this case, the RGB sample provides an OLR of 9.7 ± 0.5 kpc, which is outside the solar radius. A word of caution: because we estimate Ω 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 angular velocity of OB stars, which are younger and less affected by asymmetric drift, is slightly higher than that of RGB stars. By extrapolation to the inner disc, our Ω bar value could thus be biased towards lower values and should be taken as a lower limit.

Kinematics of the outer disc
In this section, we present and discuss the velocity maps for the OB and RGB sample in the outer disc (i.e. with galactocentric radius R > 5 kpc). We also 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 ϕ , V Z maps in the left panels. In the right panels we show the same maps (with the exception of the middle panel, where V ϕ has been replaced by the residual map of V ϕ , calculated as explained in Sec. 5.1), but compared to the spiral arms found in the overdensity (see Sec. 4), overlaid as grey shaded contours. We note that for the OB stars, the overdensity in the spiral arms and the streaming motions can be studied using the same stellar population (although the two samples are not exactly the same because only the OB stars with line-of-sight velocities were 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 ≈ 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 ϕ with respect to the mean V ϕ (R). Based on this map, the stars located just outside and inside the Local Arm move systematically more slowly than the mean V ϕ (yellow-green regions), while stars lying on the Local Arm appear to move systematically faster than V ϕ (red region). However, 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 and vertical bending waves, all with a relatively short radial wavelength.
In contrast with the OB sample, the RGB sample is an older, dynamically relaxed stellar population extending much farther from the Sun. This allows 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 masked out the region R < 5 kpc and show the V R map and the V ϕ 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 compared them with the two-arm model from Drimmel (2000), shown as solid lines, while the π/2 phase-shifted geometry, corresponding to the minimum inter-arm 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.
Moreover, the V R map (left panel) shows a large positive V R feature in the third quadrant (between 180 • < l < 270 • ) that becomes positive at about l = 170 • 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 ϕ residuals in the right panel of Figure 22 do not show any clear non-axisymmetric features in the outer disc. There is a noticeable gradient in the V ϕ residuals, however, at about y = 0, where the residuals are systematically lower for higher values of |y|. Because it is symmetric about the x−axis, this is likely a systematic in the velocity map caused by distance uncertainties, as noted in sec 3.4 (Figure 11). We also recall that features near the edge of these maps may well be artefacts from oversampling sources with overestimated 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 the (a)symmetry about the plane of the Galaxy. In Figure 23 we present velocity difference maps between the z > 0 and z < 0 hemispheres for the RGB sample (top panels) and OB stars (bottom panels). In each of the components (k = R, ϕ, Z) we created velocity maps for the upper and the lower hemispheres and then took their difference (i.e. ∆V k = V k,Z>0 − V k,Z<0 ). As with the other velocity maps, we also overplot the two-arm spiral model (black) based on 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 shows that in ∆V ϕ , 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 rotates faster below the disc plane (z < 0) by up to 10 km s −1 in the third quadrant and in the anticentre direction. A similar feature is seen in the ∆V R map, but in the opposite sense, with more positive radial motion above the disc plane. By comparing this with Figure 22 (left panel), we see that this feature corresponds with the positive V R feature in the third quadrant already noted above. Clearly, 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, ∆V R ∼ 0 km s −1 in the rest of the map. Finally, in Figure 23(c), we show the ∆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 ∆V z is about 2-4 km s −1 , that is, much lower than in the ϕ and R components. Nevertheless, the outer disc again shows a marked feature, but more towards the anticentre. We also note that at least one of the features in the ∆V z map is roughly aligned with the location of the Perseus arm from Reid et al. (2019). The ∆V map for the OB stars covers a much smaller extent of the disc, so we can only construct difference maps not much beyond 1 kpc from the disc. These are shown in the lower panels of Figure 23, with the spiral overdensity contours overlaid as in 21. In Figure 23(d), we note that the residuals in ∆V ϕ , 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 ∆V R is generally negative inside the arms, implying that the OB stars in the upper disc move towards the Galactic centre inside these arms. In Figure 21 we already noted that for the OB population, ∆V R is generally negative between the arms; the ∆V R map for the OB stars shows that this is stronger in the upper disc. Similarly, Figure 21 also showed that the ∆V ϕ residuals were also higher in the inter-arm region, and the ∆V ϕ map shows that this is more prominent in the upper disc. This suggests that there may be some shearing motion in the R and ϕ 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 ∆V z map for the OB stars in Figure 23(f). This map shows that there may be an associated breathing mode, although here we see a net expansion that may be aligned with the overdensity of the Local Arm. In general, the amplitude of the ∆V k for the OB stars is much lower than for the RGB stars.

Discussion
In section 4 we 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, P21). Segments of the nearest spiral arms are evident, although we made no attempt to derive the local stellar density. This would require an accurate 3D extinction map, as well as the selection function of the sample being used, which depends both on the survey selection function and on the criteria used to construct the sample of tracers . For instance, to compile our OB and RGB samples, we 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 and , about half of the brighter stars lack astrophysical parameters, which leads to an artificial depression in the apparent stellar density within about 1 kpc of the Sun.
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. The radial features that are so obviously visible in the distribution of our sources are only weakly visible in the velocity and dispersion maps in section 5, and they 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 overestimated distances, which in turn leads to an overestimated mean velocity component tangential to the line of sight. This causes a directiondependent distance limit to the velocity maps beyond which systematic errors dominate. These systematic errors show symmetries with respect to the line from the Sun to the Galactic centre.

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. As is the usual practice, we assumed 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 did not make any effort to distinguish between the young (thin-) disc stars and the old α-rich (thick-) disc in our RGB sample (Gaia Collaboration et al. 2022c). With these assumptions and caveats in mind, we here discuss the astrophysical significance of some of the features that we 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, we can now map a significantly larger portion of the disc based on the larger sample of stars with line-of-sight velocities. The new data reveal 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 on the azimuthal velocity field and on 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 ± 2.5 • with respect to the Sun-GC direction, as measured within 1.5 < R ≤ 2.8 kpc, and a corotation radius of 5.4 ± 0.5 kpc. Then, using the observed angular frequency, we deduce a lower limit of the bar pattern frequency of Ω bar = 38.1 +2.6 −2 km s −1 kpc −1 , and an upper limit for the location of the outer Lindblad resonance of 9.7 ± 0.5 kpc within the framework of the epicycle approximation. Bland-Hawthorn & Gerhard (2016) compiled previous estimates of the corotation radius and gave Ω bar = 43 ± 9 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±3 km s −1 kpc −1 and a corotation radius of 5.5 ± 0.4 kpc. Other works found 37.5 km s −1 kpc −1 (Clarke et al. 2019), 39 ± 3.5 km s −1 kpc −1 (Portail et al. 2017a) and 37.5 − 40 km s −1 kpc −1 (Li et al. 2022), which is based on hydrodynamical simulations and made-tomeasure models. Our findings are thus in good agreement with other estimates of the fundamental parameters of the bar. For the bar orientation angle, Bland-Hawthorn & Gerhard (2016) compiled a range of 28 • − 33 • , while other recent estimates range from 20 • to 28 • (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 disc has a large variety of features and streaming motions, which likely contain fundamental clues for 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 = 11kpc. Faure et al. (2014) carried out testparticle simulations to predict the global stellar response to spiral perturbations in the Galactic disc in the absence of external excitation (e.g. from an accreting satellite). They integrated stellar orbits in a two-arm Lin-Shu -type spiral potential (without a bar) and produced maps of mean galactocentric radial velocity (V R ). They showed (their figure 6) that inside the spiral arm corotation radius, in the region traced by the arm, the mean V R is negative (about -7 km s −1 ), that is, stars have a bulk motion towards the Galactic centre. In the region between the arms, the stellar radial motion is positive, that is, stars have a bulk motion towards the anticentre. Outside corotation, the pattern is reversed. Trends in V R induced by the spiral arms were 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 inside the dashed line in Fig.  22 (left panel), marking the location of the minimum inter-arm density.
Other features are also seen in the ∆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 disc beyond R = 10−11 kpc in the RGB sample that does not yet have a clear explanation. This feature is likely related to an asymmetry already noted in Gaia Collaboration et al. (2018b), where the authors followed the asymmetry in the azimuthal and vertical components, showing a clear bimodality of in the (V ϕ ,V z ) plane. Stars are concentrated mainly in two clumps, one with negative V z at lower V ϕ , which is more prominent in the north, and one with positive V z at higher V ϕ , which is 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 2D maps for OB stars. We recall, however, that our sample of OB stars is expected to trace the local (out to ≈ 2 − 2.5 kpc in heliocentric distance) motions of the gas, rather than the large-scale features of the Galaxy, as they inherit the motion of the gas from which they were recently born. However, because the region sampled by our OB sample is relatively small, the streaming motions in the OB stars really map only the Local Arm.
A comparison between the OB and RGB star velocity fields revealed numerous differences. The observed differences are not unexpected because these two samples trace 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.

Caveats and shortcomings
No method or procedure is perfect,and here we have adopted one that is certainly not without imperfections. In this section we confess all of our sins.
In section 5 we constructed maps of the mean velocity field in galactocentric cylindrical coordinates based on the derived velocities in the same coordinates as the individual sources. These in turn were derived from the individual measures of the proper motions and estimated distances. While our method is relatively straightforward, it has several shortcomings that should be addressed in the future. Beyond a distance of about 3 kpc, the largest source of uncertainty in our velocity maps is derived from the distance uncertainties, yet our treatment of these is not completely satisfactory.
Our distances are taken from CBJ2021, who used a Bayesian approach to estimate the distances from astrometry and photometry, incorporating a prior that includes both current astrophysical knowledge of stellar structure and some informed assumptions about Galactic structure. The probability distribution of the individual distances is asymmetric in general, but we 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 Sect. 3.3, this was done in part out of necessity as we are not able to sample the probability distribution function of the distance for each star.
We assumed no correlation between the distance and proper motion uncertainties. This cannot be true because these distances are partly 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 accounted for the correlations introduced by the coordinate transforms and distance uncertainties, which would result in a correlation between the two velocity components in galactocentric coordinates. However, we did not take these correlations into account when we estimated the mean velocity field in section 5.
Most grievously, when estimating the mean velocity and velocity dispersion for a given volume element, we only considered stars found in this small volume element, according to their estimated distances. That is, although we took the uncertainties of the distances into account when we derived the uncertainties in the individual velocities, we implicitly assumed that these distances were perfect when we binned the stars into cells to estimate the mean velocity field and velocity dispersions.
Finally, we recall 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. In addition to the CBJ2021 photogeometric distances, we attempted to use different distance estimates: the 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 an incorrect 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.

Conclusions
Using Gaia DR3 we have mapped the kinematics of the stars over an extensive area of the Milky Way disc. This has been made possible by the new line-of-sight (radial) velocities, and by new and reliable astrophysical parameters that have allowed us to distinguish 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 those of 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 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 NIR. With regard to the non-axisymmetric structure in 3D configuration space, we made no attempt to reconstruct the stellar density. Nevertheless, we were able to map the young OB sample and young clusters out to 4 to 5 kpc 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 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 regard to Galactic structure. Our kinematic maps are derived from velocities of individual stars based on individual distance estimates. A better 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. With regard 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, and the relevant parameters of the model (e.g. the bar orientation) should then be adjusted to arrive at the most likely set of parameter values given the data. Alternatively, a suite of simulations of the kinematics of the Milky Way could be generated, assuming different bar parameters, and the uncertainties modelled to transform them into a suite of mock catalogues. The one that best matches the Gaia data could then be determined. In either case, any fitting should ideally be made against the measurements in data-space rather than in model-space, and take the selection effects on the data into account.
Although we have only taken a first look at the treasure trove 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. We can certainly expect further discoveries and insights from the community in the future as it digests this latest census of the stars in our Milky Way.
In this specific case, there now is an analytic solution for σ * k . The Hessian from Eq. (B.7) reduces to (B.10) Given V k = 1 N * N * j v k, j , we have N * j (v k, j − V k ) = 0 such that the off-diagonal elements of H disappear and we are left with (B.11) Since this Hessian has a 2×2 format, we can easily invert it by hand to obtain the covariance matrix C. Indeed, we are only interested in the first element of the diagonal of C, Consequently, the uncertainty on the mean velocity V k is approximately if the dispersion satisfies σ * k ≫ σ v k, j . covered by the RGB sample. We therefore decided to fit V R as well with both mock and real kinematics. Fig. C.1 shows that the angle ϕ b that the main axis of the bisymmetry creates with Y = 0, the pseudo direction Sun − Galactic centre, can be measured independently using V R or V ϕ . This angle is that of the bar at low radius. With V R , this angle is the direction in which the right term of Eq. 20 changes its sign, that is, ϕ b = ϕ R − π/4. With V ϕ , it is the angle in which this term is minimum, that is, ϕ b = ϕ ϕ − π/2. An interesting feature the middle right panel of Fig C.1 (mock data with uncertainties) is that the direction of the change of sign of V R is no longer perfectly aligned with the major axis of the bar, but is closer to the Sun − Galactic centre line. We thus expect from the derivation of the simple asymmetric model that ϕ bisym is smaller than the true bar orientation. Figure C.2 shows the amplitudes and phase angles of the bisymmetric model of V R and V ϕ applied to the simulation with a bar phase angle of −20 • , where we masked bins with X > 0 to have the spatial distribution almost similar to that of the Gaia data. The bottom and top panels show the results obtained with and without uncertainties, respectively, to the mock variables. The asymmetric model applied to V ϕ recovers the bar orientation within R = 2.2 kpc well, regardless of the uncertainties, while for V R , the phase angle of the bar is only recovered when the uncertainties are not present, and is overestimated otherwise.
To understand why the bisymmetric model works better with V ϕ than with V R when uncertainties are added to the variables, we show in Fig. C.3 the phase angle of V R (blue curves) and V ϕ (green curves) obtained from a set of three simulations. The only difference between them is the relative uncertainty in distance used, namely 10%, 15%, and 20%. As expected, it shows that the larger the uncertainty in distance, the larger the difference between the recovered and the real phase angle. Interestingly, we also see that for V ϕ the difference between the estimated bar angle and the real orientation starts to be significant, but only in the case of the largest relative uncertainty in distance of 20%. The reason for the larger discrepancy in the radial case is that geometrically, an error on the heliocentric distance translates into an incorrect azimuthal angle, and a small change in azimuth affects the radial velocity more strongly than the azimuthal component.

Appendix D: Gaia funding acknowledgements
The Gaia mission and data processing have financially been supported by, in alphabetical order by country: -