Gaia Data Release 2 Mapping the Milky Way disc kinematics

Context. The second Gaia data release ( Gaia DR2) contains high-precision positions, parallaxes, and proper motions for 1.3 billion sources as well as line-of-sight velocities for 7.2 million stars brighter than G RVS = 12 mag. Both samples provide a full sky coverage. Aims. To illustrate the potential of Gaia DR2, we provide a ﬁrst look at the kinematics of the Milky Way disc, within a radius of several kiloparsecs around the Sun. Methods. We beneﬁt for the ﬁrst time from a sample of 6.4 million F-G-K stars with full 6D phase-space coordinates, precise par- allaxes ( σ (cid:36) /(cid:36) ≤ 20%), and precise Galactic cylindrical velocities (median uncertainties of 0.9-1.4 kms − 1 and 20% of the stars with uncertainties smaller than 1 kms − 1 on all three components). From this sample, we extracted a sub-sample of 3.2 million giant stars to map the velocity ﬁeld of the Galactic disc from ∼ 5 kpc to ∼ 13 kpc from the Galactic centre and up to 2 kpc above and below the plane. We also study the distribution of 0.3 million solar neighbourhood stars ( r < 200 pc), with median velocity uncertainties of 0.4 kms − 1 , in velocity space and use the full sample to examine how the over-densities evolve in more distant regions. Results. Gaia DR2 allows us to draw 3D maps of the Galactocentric median velocities and velocity dispersions with unprecedented accuracy, precision, and spatial resolution. The maps show the complexity and richness of the velocity ﬁeld of the galactic disc. We observe streaming motions in all the components of the velocities as well as patterns in the velocity dispersions. For example, we conﬁrm the previously reported negative and positive galactocentric radial velocity gradients in the inner and outer disc, respectively. Here, we see them as part of a non-axisymmetric kinematic oscillation, and we map its azimuthal and vertical behaviour. We also witness a new global arrangement of stars in the velocity plane of the solar neighbourhood and in distant regions in which stars are organised in thin substructures with the shape of circular arches that are oriented approximately along the horizontal direction in the U − V plane. Moreover, in distant regions, we see variations in the velocity substructures more clearly than ever before, in particular, variations in the velocity of the Hercules stream. Conclusions. Gaia DR2 provides the largest existing full 6D phase-space coordinates catalogue. It also vastly increases the number of available distances and transverse velocities with respect to Gaia DR1. Gaia DR2 offers a great wealth of information on the Milky Way and reveals clear non-axisymmetric kinematic signatures within the Galactic disc, for instance. It is now up to the astronomical community to explore its full potential.


Introduction
Our position in the disc of the Milky Way does not allow us to capture the global picture of our galaxy easily. Mapping its 3D structure requires large and precise astrometric catalogues. The second Gaia data release (Gaia DR2, Gaia Collaboration 2018) contains positions and parallaxes for 1.3 billions sources down to magnitude G ∼ 21 mag, which multiplies by a huge factor the number of stars for which a distance can be derived with respect to Gaia DR1. Not only does Gaia DR2 provide the 3D location of a very large sample of stars in the Galaxy, it also contains full velocity information (proper motions and line-of-sight velocity) for 7.2 million stars brighter than G RVS = 12 mag, and transverse velocity for an unprecedentedly large number of stars. This paper belongs to a series of six Gaia DR2 performance verification papers that are meant to demonstrate the quality of the catalogue through a basic examination of some of the key science cases of the Gaia mission. In this paper, we report a first look at the kinematic properties of the Milky Way disc as pictured by the second Gaia data release.
Non-axisymmetric features manifest themselves not only in configuration spaces, but also in kinematic spaces, where they leave specific signatures related to their spatial extension, rotation speed around the Galaxy centre, and growth rate (Siebert et al. 2012;Faure et al. 2014;Monari et al. 2014Monari et al. , 2016bDebattista 2014;Bovy et al. 2015;Grand et al. 2015Grand et al. , 2016Antoja et al. 2016;Pasetto et al. 2016). Many studies prior to Gaia (Eggen 1958(Eggen , 1996Chereul et al. 1999;Dehnen 1998;Famaey et al. 2005; Antoja et al. 2008; Gómez et al. 2012a), and especially since the epoch of the HIPPARCOS satellite (Perryman et al. 1997), have studied the kinematics of stars in the solar neighborhood and have shown that the stellar velocity and phase-space distributions are not smooth, but rather clumpy. Several hypotheses were able to explain the nature of this clumpy distribution, suggesting that they might be remnants of stellar clusters (Eggen 1996), substructures related to orbital effects of the bar and/or the spiral arms (e.g. Dehnen 2000;De Simone et al. 2004;Quillen & Minchev 2005;Chakrabarty 2007; Antoja et al. 2009), remnants of accreted systems from the halo (Helmi et al. 1999;Villalobos & Helmi 2009;Gómez & Helmi 2010;Re Fiorentin et al. 2015;Jean-Baptiste et al. 2017), or substructures induced in the stellar disc by external perturbations Minchev et al. 2009;Gómez et al. 2012b;Jean-Baptiste et al. 2017). Despite all this theoretical and observational work, it is still an open issue how we can distinguish between the different types of substructures. With RAVE (Steinmetz et al. 2006), LAMOST (Liu et al. 2014) combined with TGAS (Gaia Collaboration 2016; Lindegren et al. 2016), and APOGEE-2 South (Majewski et al. 2016(Majewski et al. , 2017, Antoja et al. (2012Antoja et al. ( , 2014, Monari et al. (2017) and Hunt et al. (2018) concluded that at least one of these substructures, the Hercules stream, evolves with Galactic radius, consistently with the effects of the Outer Lindlblad Resonance of the bar. However, other studies have suggested a pattern speed for the Milky Way bar that is slower than previous estimates, placing this resonance well outside the solar radius (Liu et al. 2012;Portail et al. 2017;Pérez-Villegas et al. 2017). To understand the role of the stellar bar, it is necessary both to map the kinematics of disc stars in the Galaxy over a larger spatial extent and to increase the statistics (the number of stars with full 3D kinematic information) out to a few kpc from the Sun. Extending the spatial scale of kinematic studies to larger regions of the Galactic disc is also essential for quantifying the amplitude of velocity gradients, detection of which is now limited to a region of a few kiloparsec around the Sun (see Siebert et al. 2011;Carrillo et al. 2018;Liu et al. 2017), and constrain their origin.
In addition to secular evolutionary processes, a disc galaxy like ours is expected to have experienced several accretion events in its recent and early past (Bullock & Johnston 2005;De Lucia & Helmi 2008;Stewart et al. 2008;Cooper et al. 2010;Font et al. 2011;Brook et al. 2012;Martig et al. 2012;Pillepich et al. 2015;Deason et al. 2016;Rodriguez-Gomez et al. 2016). While some of these accretions are currently being caught in the act, like for the Sagittarius dwarf galaxy (Ibata et al. 1994) and the Magellanic Clouds (Mathewson et al. 1974;Nidever et al. 2010;D'Onghia & Fox 2016), we need to find the vestiges of ancient accretion events to understand the evolution of our Galaxy and how its mass growth has proceeded over time. Events that took place in the far past are expected to have induced a thickening of the early Galactic disc, first by increasing the inplane and vertical velocity dispersion of stars (Quinn et al. 1993;Walker et al. 1996;Villalobos & Helmi 2008Zolotov et al. 2009;Purcell et al. 2010;Di Matteo et al. 2011;Qu et al. 2011;Font et al. 2011;McCarthy et al. 2012;Cooper et al. 2015;Welker et al. 2017), and second by agitating the gaseous disc from which new stars are born, generating early stellar populations with higher initial velocity dispersions than those currently being formed (Brook et al. 2004(Brook et al. , 2007Forbes et al. 2012;Bird et al. 2013;Stinson et al. 2013). These complementary modes of formation of the Galactic disc can be imprinted on kinematics-age and kinematics-abundance relations (Strömberg 1946;Spitzer & Schwarzschild 1951;Nordström et al. 2004;Seabroke & Gilmore 2007;Holmberg et al. 2007Holmberg et al. , 2009Bovy et al. 2012aBovy et al. , 2016Haywood et al. 2013;Sharma et al. 2014;Martig et al. 2016;Ness et al. 2016;Mackereth et al. 2017;Robin et al. 2017), and distinguishing between them requires full 3D kinematic information for several million stars, in order to be able to separate the contribution of accreted from in-situ populations, and to constrain impulsive signatures that are typical of accretions ) versus a more quiescent cooling of the Galactic disc over time. Accretion events that took place in the more recent past of our Galaxy can also generate ripples and rings in a galactic disc (Gómez et al. 2012b), as well as in the inner stellar halo (Jean-Baptiste et al. 2017). Such vertical perturbations of the disc are further complicated by the effect of spiral arms Monari et al. 2016b), which together with the effect of accretion events might explain vertical wave modes as observed in SEGUE and RAVE (Widrow et al. 2012;Williams et al. 2013;Carrillo et al. 2018), as well as in-plane velocity anisotropy (Siebert et al. 2012;Monari et al. 2016b). Mapping the kinematics out to several kiloparsec from the Sun is crucial for understanding whether signs of these recent and ongoing accretion events are visible in the Galactic disc, to ultimately understand to what extent the Galaxy can be represented as a system in dynamical equilibrium (Häfner et al. 2000;Dehnen & Binney 1998), at least in its inner regions, or to recover the nature of the perturber and the time of its accretion instead from the characteristics and strength of these ringing modes (Gómez et al. 2012b).
Signatures of interactions and gravitational disturbances of satellite galaxies can also affect the outer disc beyond the solar radius, in regions where the stellar surface density drops and the disc is more fragile to external perturbations. Several works have discussed the possibility that the Galactic warp may be generated by the interaction with the Magellanic Clouds (Burke 1957;Weinberg & Blitz 2006) or Sagittarius (Bailin 2003), while other scenarios suggest that a warped structure in a galaxy disc may be generated by a dark matter halo distribution that is off-centred or tilted with respect to the baryonic one (Bailin & Steinmetz 2003), by bending instabilities (Revaz & Pfenniger 2004) in the disc, or by misaligned infall of material (Ostriker & Binney 1989;Quinn & Binney 1992). These scenarios predict either long-lived, transient, or repeatedly excited structures, and it is clear that to understand the origin of the Galactic warp, we need to understand its dynamical nature, since, for example, a long-lived warp would leave a specific signature in the kinematics of stars in the outer disc (Abedi et al. 2014;Poggio et al. 2017).
In the coming years, the astronomical community will work towards answering these great questions about the Galaxy with the help of Gaia data. In this paper, we provide a first exploration of the kinematic properties of the Milky Way disc that already reveals novel results, shows the far-reaching possibilities of the data, and predicts their high future impact. The paper starts by a description of the Gaia DR2 data that are used in this analysis (Sect. 2). Details are given about calculating distances, velocities, and their uncertainties, as well as about the different data selections. In Sect. 3 we start by exploring the velocity components in 3D, their medians and dispersions, by searching for global trends as a function of position, distance from the Galactic centre, and height above the plane. This analysis for the first time presents full 3D kinematic maps of the Galaxy up to several kiloparsec from the Sun. In Sect. 4 we zoom into the solar neighborhood and revisit its velocity distribution by searching for kinematic substructures at small scales with unprecedented accuracy, and also by showing how they evolve with spatial position. The full-sky coverage of Gaia overcomes limitations in angular coverage of earlier studies. Finally, in Sect. 5, we present the main conclusions of this work.

Data
In this section, we describe and characterise briefly the Gaia DR2 data that we used. We start with an overview of the content of DR2. Secondly, we detail how the distances, velocities, and their uncertainties are calculated. Next, we explain how we built a dereddened HR diagram to select different stellar populations, followed by details on the different data samples that are used throughout the paper, and details on their main characteristics. Finally, the last two subsections briefly characterise important aspects of the samples, such as the correlations between variables, and the anisotropy of the samples.

DR2 data overview
Gaia DR2 provides astrometric parameters (positions, parallaxes, and proper motions) for 1.3 billion sources. The median uncertainty for the bright sources (G < 14 mag) is 0.03 mas for the parallax and 0.07 mas yr −1 for the proper motions. The reference frame is aligned with the International Celestial Reference System (ICRS) and non-rotating with respect to the quasars to within 0.1 mas yr −1 . The systematics are below 0.1 mas and the parallax zeropoint uncertainty is small, about 0.03 mas. Significant spatial correlations between the astrometric parameters are also observed. For more details about the astrometric content of Gaia DR2, see Lindegren et al. (2018), Arenou et al. (2018) and references therein.
The photometric content of Gaia DR2 consists of weightedmean fluxes and their uncertainties for three passbands, G, G BP , and G RP . All sources have G photometry, but only about 1.4 out of the 1.7 billion sources have both G BP and G RP photometry. The sources without colour information mainly lie in crowded regions where the larger windows for the BP and RP photometers have a higher chance of overlap between sources and make the photometry unreliable. The processing for future data releases will include deblending algorithms that will increase the number of sources with colour information. The precision at G = 12, the magnitude most relevant for this kinematic study, is around 1 mmag or better for all three passbands. However, there are systematics in the data at the 10 mmag level. For more details about the photometric content of Gaia DR2, see Evans et al. (2018) and references therein.
To facilitate the selection of specific types of stars, we also used the extinction A G and color excess E(G BP − G RP ) provided in Gaia DR2, whose estimation was described in Andrae et al. (2018). However, the accuracy of the astrophysical parameters, derived from Gaia data alone, is degenerate for some parts of the Hertzsprung-Russel (HR) diagram, especially for high extinction values. To assist in the sample selections, we therefore also made use of 2MASS photometry of the Gaia sources, specifically, of the Gaia/2MASS cross-match provided within GACS for Gaia DR2 (see Marrese et al. 2018). Details of how the 2MASS photometry was used are described below in Sect. 2.3 and 2.4.
A novelty of Gaia DR2 with respect to Gaia DR1 is that it contains line-of-sight velocities 1 for 7.2 million stars brighter than G RVS = 12 mag that were observed with the Radial Velocity Spectrometer . The stars are distributed throughout the full celestial sphere. This release contains line-ofsight velocities for stars with effective temperatures in the range ∼[3550, 6900] K. Cooler and hotter stars will be published in future Gaia releases. The precision of Gaia DR2 line-of-sight velocities is at the km s −1 level. At the bright end, the precision is of the order of 0.2 to 0.3 km s −1 . At the faint end, it is of the order of 1.4 km s −1 for T eff = 5000 K stars and ∼3.7 km s −1 at T eff = 6500 K. For more details about the Gaia spectroscopic processing pipeline and the Gaia DR2 line-of-sight velocities, see Sartoretti et al. (2018) and Katz et al. (2018) and references therein.
The global validation of Gaia DR2 is described in Arenou et al. (2018) and references therein.

Calculation of distances, velocities, and uncertainties
In order to map the stars in position and velocity space, we must derive distances from the Gaia astrometry. For this purpose, we have selected only stars with /ε > 5 and adopted 1/ as our distance estimate. It is well-known that the inverse of the parallax is biased when the uncertainty in parallax is significant (Brown et al. 1997;Arenou & Luri 1999;Luri et al. 2018). To quantify the distance bias introduced when using 1/ as a distance estimator and a cut at 20% relative uncertainty in parallax, we used the simulations described in Sect. 2.4. We established that inverting the parallax leads to unbiased distances out to about 1.5 kpc, with overestimates of the order of 17% at 3 kpc. We therefore have to bear in mind that the distance bias in the extremes of our main sample is non-negligible.
Note that this cut in relative uncertainty in parallax results in a cut in apparent magnitude, and other minor selection effects might be caused by this. However, after tests with our set of simulations, we concluded that this cut does not introduce relevant artefacts in the kinematics. Alternatively, Bayesian methods might be used to infer distances from parallaxes instead of selecting stars with small relative uncertainty (e.g. Bailer-Jones 2015). However, this is more complex in the sense that they require fixing a prior, and even the simplest sensible prior involves numerical solutions for most estimators and for all the confidence intervals. In this exploratory study, we chose to select small uncertainty in parallax since it is simpler and serves the purposes of our work well.
Gaia provides the five-parameter astrometric solution 2 and line-of-sight velocities, (α, δ, , µ * α , µ δ , V los ), together with their associated uncertainties and correlations between the astrometric quantities. From these observables and the derived distances, we computed heliocentric and Galactic Cartesian and cylindrical positions and velocities. For the Cartesian heliocentric velocities, we took the usual convention of U, V, and W oriented towards the Galactic centre, the direction of Galactic rotation, and the north Galactic pole, respectively. The Galactic cylindrical coordinates are (R, φ, Z, V R , V φ , V Z ) with φ in the direction of Galactic rotation and with an origin at the line Sun-Galactic centre. The Cartesian Galactic coordinates are oriented such that the Sun is located at the X negative axis. For these transformations, we needed to adopt a height of the Sun above the plane. We used the value given by Chen et al. (2001) of 27 pc, although other values can be 14 ± 4 pc from 2 Proper motion in right ascension µ * α ≡ µ α cos δ of the source in ICRS at the reference epoch. This is the projection of the proper motion vector in the direction of increasing right ascension. COBE/DIRBE (Binney et al. 1997) or 15.3 +2.24 −2. 16 from Gaia DR1 (Widmark & Monari 2017). We also adopted the distance of the Sun to the Galactic centre R of 8.34 kpc and the circular velocity at the solar radius of V c = 240 km s −1 from Reid et al. (2014). We took the peculiar velocity of the Sun with respect of the local standard of rest from Schönrich et al. (2010), that is, (U , V , W ) = (11.1, 12.24, 7.25) km s −1 . The resulting value of (V c + V )/R is 30.2 km s −1 kpc −1 , which is compatible with the value from the reflex motion of Sgr A* of Reid & Brunthaler (2004). In these coordinate transformations, we propagated the full covariance matrix. This means that we have the correlations between uncertainties in Cartesian and cylindrical coordinates at our disposal.

Intrinsic colour computation
To select stars in the HR diagram, we have used cuts in absolute magnitude and intrinsic colours. For this an extinction correction needed to be applied, in particular for distant giants and hot stars. While first extinction estimates by Gaia consortium have been made using the Gaia integrated bands alone, the addition of the 2MASS colours strongly helps to break the T eff -extinction degeneracy (Andrae et al. 2018). We used here the Gaia DR2 provided cross-match with 2MASS (Marrese et al. 2018). We selected 2MASS stars with photometric quality flag AAA and photometric uncertainties lower than 0.05 mag. We used the same Gaia photometric cuts as in Babusiaux et al. (2018): photometric uncertainties smaller than 5% for G BP and G RP and 2% for G, and a selection on the G BP /G RP excess flux factor based on the star colour. To derive intrinsic colour-colour relations, we selected low-extinction intrinsically bright stars as in Babusiaux et al. (2018), for example, using the 3D extinction map of Capitanio et al. (2017) 3 and the Gaia DR2 distances, to select stars with E(B − V) < 0.015 and M G < 2.5. For each photometric band X = G BP , G RP , J, H, we built a fifth-order polynomial relation to model (G − X) 0 as a function of (G − K s ) 0 . We used the extinction coefficient models described in Danielski et al. (2018), computed using the nominal passbands. We pre-selected intrinsically bright stars using the 2MASS K s magnitude, which is less strongly affected by extinction: where the astrometry is given in milliarseconds. Then the extinction A 0 and the intrinsic colour (G − K s ) 0 were determined for each star through a maximum likelihood estimator (MLE). This takes into account the photometric uncertainties, the intrinsic scatter around the intrinsic colour-colour relation (which is between 0.01 and 0.03 mag), and the validity intervals of these relations as well as the positivity of the estimated extinction. A chi-square test was performed to verify the validity of the resulting parameters, removing stars with a p-value limit lower than 0.05. We also removed stars for which the MLE did not converge and those with an error on (G − K s ) 0 larger than 0.5 mag. In total, we obtained extinction corrections for 90% of the sample. Figure 1 shows the de-reddened HR diagram.

Data selection
As discussed above in Sect. 2.2, we selected sources with /ε > 5. This cut selects stars with positive parallaxes and a relative parallax uncertainty smaller than 20%. After this cut, we further selected several samples that we use in the different sections of this study.
1. Main sample. This sample consists of the 6 376 803 sources with an available five-parameter astrometric solution, line-of-sight velocities, and /ε > 5. The intrinsic magnitudes and colours were calculated using Gaia and 2MASS photometry, as explained in Sect. 2.3. In the top and bottom panels of Fig. 2, we show the surface density per bins of 100 pc × 100 pc in (X, Y) and (R, Z) planes, respectively, while in top and bottom panels of Fig. 3, we show the G apparent magnitude and the Galactic radius distribution of the main sample (black lines) and the remaining working samples. For these stars, we computed the full 6D phase space coordinates as detailed in Sect. 2.2. The top panel of Fig. 4 shows the distribution of uncertainties in Galactic cylindrical velocities of the main sample. The median uncertainties are (ε V r , ε V φ , ε V z ) = (1.4, 1.4, 0.9) km s −1 , and 20% of the stars have an uncertainty in all velocity components that is smaller than 1 km s −1 . The distributions in ε V r and ε V φ are similar and differ from the distribution for ε V Z , which is more precise. The reason is that most of the stars are located in the Galactic plane: for these stars, the main contribution to the vertical velocity comes from the astrometric quantities, which for this sample have smaller uncertainties than does the line-of-sight velocity. The uncertainties as a function of distance are shown in the bottom panel. They seem to increase approximately linearly in this log-log plot. The median velocity uncertainty is below 1 km s −1 at distances closer than 0.5 kpc, and below 2 km s −1 at distances closer than 2 kpc. In addition, uncertainties larger than 10 km s −1 are only reached at distances larger than 5 kpc.
The main sample supersedes any previous full 6D phasespace sample in terms of quantity and precision of the data. For instance, the main sample is about 12 times larger in number of stars than a sample made from UCAC proper motions (Zacharias et al. 2013) and RAVE line-of-sight velocities and derived spectro-photometric distances (Kunder et al. 2017). Thus, the statistics enable studying the Galaxy kinematics in more details and at much larger distances than before. At the faint end, the precision of the RAVE line-of-sight velocities is comparable to that of the RVS. However, with Gaia DR2, the precisions as a function of distance in the derived distances and in the proper motions are about two and more than ten times better, respectively. This combination means that the precision in Galactocentric cylindrical velocities of the main sample is approximately 5-7 times better. As an example, we show the Toomre diagram of the main sample in Fig. 5.
2. Giant sample. This is a sub-selection of the main sample that includes only giant stars selected on their absolute magnitude in G band M G < 3.9 and intrinsic colour (G BP − G RP ) 0 > 0.95. The intrinsic magnitudes and colours were calculated by using Gaia and 2MASS photometry, as explained in Sect. 2.3. This sample contains 3 153 160 sources. As noted in Fig. 3, about half of the stars in the main sample are (red) giants, which are the main contribution at distances larger than 1 kpc from the Sun. That is why this sample is used in Sect. 3 to analyse the large-scale kinematic maps in the Galactic disc. As expected, 78% of the (red) giant sample is located within 3 kpc of the Sun. Nonetheless, the inner regions, that is, areas towards the Galactic centre with Galactic radius between 3-5 kpc, are still well sampled with more than 500 000 stars (see the bottom panel of Fig. 3). Furthermore, in the outskirts of the galactic disc, our red giant sample contains more than of the red giant branch, and their stellar evolutionary stage is therefore different from the red clump sources, most of which are located at about ±2 kpc from the Sun. The median uncertainties are (ε V r , ε V φ , ε V z ) = (1.6, 1.7, 1.2) km s −1 , and 13% of the stars have an uncertainty in all velocity components that is smaller than 1 km s −1 .
3. Solar neighbourhood sample. This is a sub-selection of the main sample with stars located within 200 pc of the Sun, that is, with > 5 mas. This comprises 366 182 stars with a median velocity uncertainty of (ε U , ε V , ε W ) = ( 0.4, 0.4, 0.4 ) km s −1 and with 78% of stars having uncertainties smaller than 1 km s −1 in all components.
4. OB sample. This is the selection of OB stars used in Sect. 3 to map the median vertical velocity of young stellar populations. This sample is different from those described above in that it is not constrained to sources with available line-of-sight velocities. However, the additional challenge is identifying young, intrinsically blue stars near the Galactic plane that are significantly reddened.
An initial list of OB star candidates in DR2 was found using the following criteria: where A G and E(G BP − G RP ) are the extinction and colour excesses provided in Gaia DR2 (see Andrae et al. 2018), and  The traditional use of the Toomre diagram to classify stars into stellar populations is complicated by the great range of the Galactic radius of the sample (Fig. 3) and the possibility that both the mean V φ of the thin disc and the V φ lag between the thin and thick disc may vary with Galactic radius. Nevertheless, it shows that the sample is dominated by the thin disc. In the solar neighbourhood, the thin disc has an azimuthal velocity close to the LSR, and the thick disc lags behind by a few tens of km s −1 .
A11, page 7 of 40 A&A 616, A11 (2018) is expressed in mas. To ensure that our sample indeed consists of young stars rather than giants or red clump stars with erroneous extinctions, a further selection was made using the 2MASS photometry that also satisfies the following conditions: These colour-colour selection criteria were adopted from those described by Poggio et al (in prep.) and are based on the observed 2MASS colours of spectroscopically bona fide OB stars from the Tycho-2 stars found in Gaia DR1 and the Tycho-2 spectral type catalogue (Wright et al. 2003). In addition, the photometric quality conditions ε J,H,K s < 0.05 and 2MASS photometric flag equal to AAA were applied to avoid sources with problematic photometry. These selections yielded 285 699 stars whose 2MASS/Gaia colours and astrometry are consistent with our sources being OB stars. However, given the relatively large uncertainties on the individual extinction parameters, our sample is likely to also contain a significant number of upper main-sequence A stars. Nevertheless, such stars, being young, still serve our purpose here. The apparent magnitude and galactocentric radial distribution is shown in Fig. 3.

Simulation of Red Clump disc stars
In order to analyse the effect of errors and biases throughout the different sections of this study, we used the simulation of Gaia data provided in Romero- Gómez et al. (2015). This is a test-particle simulation of Red Clump disc stars that evolved in a barred galactic potential. We only kept stars with G ≤ 13 from the entire simulation to mimic the magnitude distributions of our main sample. This led to a simulation of one million Red Clump disc stars with astrometric and line-of-sight velocity uncertainties that matched those of Gaia DR2. We rescaled the end-of-mission astrometric uncertainty prescribed on the Gaia Science Performance webpage (see also de Bruijne et al. 2014) to the Gaia DR2 uncertainty for 22 months of mission 4 , and for the bright stars, we included a multiplying factor of 3.6 to match the distribution of the uncertainty as a function of G magnitude observed in the Gaia DR2 data. The line-of-sight velocity uncertainties were also rescaled to match the uncertainty for the Red Clump-type of stars observed in our Gaia sample.

Correlations between astrometric and derived quantities
In Fig. 6 we show for the main sample the correlation coefficient between the Galactic radius and the different components 4 http://www.rssd.esa.int/doc_fetch.php?id=359232 of the Galactic velocity as a function of the Galactic longitude. Most of the stars are concentrated in regions of correlations near unity, which are positive or negative depending on the Galactic longitude. This behaviour is mainly due to a geometric effect and not to especially strong correlations between the Gaia observables. The stars with correlation coefficients near to 1 in these panels do not have strong correlations between the Gaia observables. We note that the median absolute correlations of this sample are ρ π−µ α = −0.03, ρ π−µ delta = 0.01 and ρ µ α −µ δ = 0.01, and for 89% of the stars, all three correlations are weaker than 0.4. The behaviour in these panels arises because both the Galactic radius and the velocities are dependent on the heliocentric distance, which in our study we take as the inverted parallax. In this sense, any uncertainty in distance would translate into a proportional uncertainty in R and (V R , V φ , V Z ), its sign depending on the position in the Galaxy. Therefore, the uncertainties in radius and velocities are highly correlated.
While the correlations on the observables might bias some derived quantities, this will only happen in the limit of large uncertainties and depending on the problem under study. We also note that if the errors on the astrometric basic parameters are random, as expected, these high correlations do not necessarily translate into a bias, meaning that this is not equivalent to having a systematic error. However, we emphasize that correlations are important in the uncertainty propagation and should not be neglected.
In our data selection we did not perform any cut in velocity uncertainty. Figure 7 shows the uncertainty in velocity as a function of velocity for the three Galactic components. Since the velocities and their uncertainties are correlated, removing stars with large uncertainties, such as those above the dashed black line at 2 km s −1 , entails the removal of the stars with higher velocities. This can cause large biases on derived quantities such as the velocity dispersion, and we have checked that even the mean velocities as a function of Galactic radius or height above the plane appear to be highly biased (with differences of up to 20 km s −1 ) when performing these data selections (see Appendix B).

Magnitude limit and asymmetric extinction
Even though Gaia is unique in covering the whole sky, the effects of the scanning law, extinction, and other complex aspects of the completeness of the data (see Arenou et al. 2018 andKatz et al. 2018) complicate the selection function. As a consequence, the properties of the main sample depend strongly on the direction. To show one example, the average vertical position Z in the X-Y plane of the giant sample is displayed Fig. 7. Median of the uncertainty in the three Galactic cylindrical velocity components (radial V R , azimuthal V φ − V c , and vertical V z ) as a function of the corresponding velocity components for the main sample. The colour-shaded areas show the 25% and 75% quartiles. The horizontal dashed black line indicates the bias that would be introduced if a cut of 2 km s −1 were performed.
in Fig. 8 (top). The median vertical position is a strong function of Galactic longitude, which is clearly affected by the extinction in our Galaxy, which is highly non-uniform. The values of median Z are higher than 600 pc at distances beyond 3 kpc. In the bottom panel of Fig. 8, the same quantity is shown for the simulation of Red Clump stars described above. In this simulation, the 3D extinction model of Drimmel & Spergel (2001) was used. Similar trends are shown between the Gaia data and the simulation. To reduce the bias on the median Z as a function of Galactic radius significantly, in Sect. 3, we divide the disc into layers of 400 pc height when it is observed face-on.
On similar lines, the uncertainties on the derived quantities also depend strongly on the position in the Galaxy in a complex way that is greatly related to extinction. Figure 9 shows the median velocity uncertainties as a function of position in configuration space. While the uncertainties globally increase as a function of distance from the Sun, as expected, this increase depends on the direction because it is affected by interstellar extinction. For instance, some blue spikes appear in these panels in lines of sight with lower extinction, while in other directions, the uncertainty achieves high values at close distances to the Sun. However, we note that the median velocity uncertainties are very small compared to other previous catalogues: they are of the order of 6-10 km s −1 only at the extremes of the sample. We also emphasise that given the large number of stars, the uncertainties on the median velocities in a given Galactic position are much smaller than these median (individual) velocity uncertainties showed here. For instance, median velocities at 1 and 1.5 kpc have unprecedented precisions of 0.5 and 1 km s −1 , respectively (see colour-shaded areas in Figs. 12-14).

Mapping the disc median velocities and velocity dispersions
Non-axisymmetric structures (e.g. bar and spiral arms) and external perturbers (e.g. the Sagittarius dwarf galaxy, the Magellanic Clouds, and dark matter sub-halos) are expected to disturb the Milky Way velocity field. In the past decade and thanks to large spectroscopic surveys and proper motion catalogues, RAVE ( Schönrich & Dehnen 2018). In this section, we take advantage of the large data volume, full sky coverage, accuracy, and precision of Gaia DR2 to re-examine these kinematic features at higher accuracy than ever before. We study the kinematics of the sample of giant stars (described in Sect. 2.4), and map the medians (Ṽ R ,Ṽ φ ,Ṽ Z ) and the dispersions (σ V R , σ V φ , σ V Z ) of the Galactocentric velocities as a function of the location in the Galaxy (X, Y, R, φ, Z).

Method
Four projections were used to study the kinematics (median velocities and velocity dispersions) of the giant sample. The layers were centred on the Galactic mid-plane, except when we studied the median vertical velocity, for which specific attention was given to the possible northsouth asymmetries. In this specific case, the giant sample was split into six layers, three above and three below the mid-plane. Each sub-sample was then divided into R-cells of 400 pc. -Galactocentric cylindrical vertical projections. The sample was first split into four azimuth slices of 15 degrees each and into three ranges in Galactic radius: [6,8], [8,10], and [10, 12] kpc. Each sub-sample was then divided into Z-cells of 200 pc. This projection was used only to study the median vertical velocity,Ṽ Z . When the cells were sufficiently populated, the medians (Ṽ i , i ∈ {R, φ, Z}) and the dispersions (σ V i , i ∈ {R, φ, Z}) of the velocities and their associated uncertainties were derived 5 . A minimum of 30 stars per cell was required to compute the moments of the velocities in the XY-maps and RZ-maps. The minimum was 50 stars for the radial projections. Each face-on or edge-on map had its own colour range dynamics in order to heighten the contrast between the spatial structures within the map. Conversely, the different layers and azimuth slices shared the same scale in the R-projections in order to facilitate the comparison. 5 According to Formulae A.1 to A.5 (see Appendix A).
The maps are (roughly) centred on the Sun (X, Y) or (R, Z) position, and the Galactic centre is located on the left side. In the face-on maps, the Milky Way rotates clockwise. Figures 10 and 11 present the face-on and edge-on views of the median velocities and velocity dispersions for the mid-plane layer. For clarity, the full mosaics of face-on and edge-on maps, which offer vertical and azimuthal tomographic views of the disc kinematics, are presented in Appendix C.
To quantify and visualise the respective contributions of bending and breathing modes, we also map the bending and breathing velocities (Fig. C.6). We calculated them as the halfsum (mean) and half-difference of the median vertical velocities in symmetric layers with respect to the Galactic mid-plane: and whereṼ Z ((X, Y), L) is the median vertical velocity in the cell (X, Y) and in the horizontal layer L. Layer L was chosen to lie in the north Galactic hemisphere, and layer −L is the symmetric layer in the south Galactic hemisphere. Formulae 7 and 8 are similar to those defined by Widrow et al. (2014), except that we calculated the half-difference for the breathing velocity, while they used the full difference. Figure 12 shows the median radial velocity,Ṽ R , as a function of Galactic radius for negative (left) and positive (right) azimuths and for different Z layers (the different curves). The median radial velocity has a U-shape, with a minimum at about 9 kpc. Around this minimum and within a broad layer below and above the mid-plane, the median radial velocity is negative, meaning that more stars move inwards than outwards. At a distance from the minimum of 1 to 2 kpc, the median radial velocity becomes positive, meaning that more stars move outwards than inwards. At negative azimuths, the median radial velocity reaches maxima at around 6.5-7.5 kpc and 11-13 kpc and then decreases again. More than a U-shape, at negative azimuths, the  median radial velocity seems to present oscillations. At positive azimuths, the signal is partially different. A maximum may be indicated at around 12-13 kpc for Z in [−1000, −200] pc, while in the other layers, the median radial velocity seems to continue to increase with Galactic radius, but the data there are too noisy to conclude. In the inward direction, the radial velocity shows a plateau starting at around 7 kpc. Farther inward, the rising of the green and orange curves at R ≤ 5−6 kpc should be considered with caution, as it is significant only at the ∼1σ level. The vertical behaviour of the radial velocity oscillation varies with azimuth. At negative azimuths, the median radial velocity is minimum in the Galactic mid-plane and increases with distance to the plane. At positive azimuths, the median radial velocity shows a much smaller vertical gradient.

Radial velocity
Seen face-on (upper left panel of Fig. 10 and Fig. C.1), the negative radial velocities (blue-green pattern) have a semicircular geometry with a small pitch angle that does not seem to present vertical variations.
Using RAVE data, Siebert et al. (2011) measured a negative radial velocity gradient from about 2 kpc inward of the Sun to about 1 kpc outward. This gradient was confirmed and further studied by several groups (Williams et al. 2013;Carrillo et al. 2018). Carrillo et al. (2018) also observed the onset of a positive gradient beyond the solar radius. Using samples of LAMOST giants, Tian et al. (2017) and Liu et al. (2017) also measured positive radial velocity and line-of-sight velocity gradients in the direction of the galactic anti-centre, which flatten at around 2 kpc beyond the Sun. Carlin et al. (2013Carlin et al. ( , 2014 studied the motions of F-type stars observed with LAMOST in the direction of the Galactic anti-centre. They observed an inward mean motion of the stars in the Galactic plane and an inversion of the sense of the mean motion at a distance from the plane, in particular at Z −0.8 kpc. The negative and positive gradients revealed by previous studies are well visible in Gaia DR2 data as part of oscillation(s) on a kiloparsec scale. The full-sky coverage and large statistics of the Gaia DR2 catalogue allows us to map the oscillation in 3D and to observe its semi-circular geometry, with a small pitch angle. At negative azimuth and around R = 9 kpc, the sign of the median radial velocity changes, that is, it is negative for |Z| 0.6−0.8 kpc and positive at larger distances from the plane (see  noted that the vertical variation ofṼ R is relatively modest, of the order of 5-10 km s −1 . Therefore a small change in the radial velocity zeropoint and in particular in the peculiar radial velocity of the Sun can modify the position of the inversion of the radial mean motion. Different methods can indeed lead to estimates of the solar peculiar radial velocity with respect to the LSR that differ by a few km s −1 : that is, U = 11.1 km s −1 (Schönrich et al. 2010) and U = 14.0 km s −1 (Schönrich 2012). Figure 13 shows the Milky Way stellar median rotation profiles from 4 to 13.2 kpc from the Galactic centre. In the inner part of the Galaxy, the median azimuthal velocity presents a steep positive gradient with Galactic radius before it reaches a maximum at around 230 km s −1 (a few km s −1 below 6 the value adopted in this study for the LSR: i.e. 240 km s −1 ). When the maximum is reached, the azimuthal velocity presents a relatively flat profile, with variations of a few km s −1 with Galactic radius. The asymmetric drift is expected to play a major role in the increase of the median velocity for increasing radius. At inner radii, the velocity dispersion in the radial velocity is larger (see Sect 3.5 and Fig 16), and the asymmetric drift correction is proportional to this dispersion squared. A detailed correction for the asymmetric drift is beyond the scope of this study, but only when this is completed can we assess whether the gradient in the azimuthal velocity is related to a gradient in the potential, to the effects of the non-axisymmetric perturbations such as the Galactic bar, and/or to the increasing weight of the thin disc with respect to the α-element-rich thick disc (the former presenting a greater radial scale length and a faster rotation than the latter, see Bovy et al. 2012b;Robin et al. 2014).

Azimuthal velocity
The rotation profiles reach their maximum at shorter radius in the mid-plane than at larger distances from the plane: R ∼ 6−7 kpc for Z in [−200, 200] pc, R ∼ 8 kpc for |Z| in [200,600] pc, and R ∼ 9−11 kpc for |Z| in [600, 1000] pc. The comparison of the two panels of Fig. 13 and the comparison of the red and orange curves, on the one hand, and of blue and purple curves, on the other hand, show that the rotation profiles are relatively symmetric in azimuth and with respect to the midplane. At R = 12 kpc, most curves are contained within a narrow range of median V φ . The decrease with radius of the vertical gradient in azimuthal velocity is also visible in the edge-on maps ( Fig. 11 middle left panel and Fig. C.4) as an outward flaring of the iso-velocity contours. This can be explained by an increase  Fig. 12 for the median vertical velocity,Ṽ Z . The disc has been divided into six layers (the six curves), three above and three below the mid-plane.
in asymmetric drift with Z. This change can be due to the different relative proportion of the thick and thin disc and/or of the different mean populations (young versus old), and to the variation in radial force in the galactic disc. Bienaymé et al. (2015) have developed a dynamically self-consistent Staeckel potential using the mass distribution of the Besançon Galaxy model. They showed the variation in asymmetric drift as a function of Galactocentric radius and distance to the Galactic plane. While the variations with R are mild between 5 and 10 kpc (less than 20%), the effect in Z is very noticeable for the thin and thick discs both. These variations are shown in Robin et al. (2017) and compared with the kinematics in Gaia DR1. The lag, typically of 5 to 20 km s −1 in the Galactic plane, can be increased by 50 to 100% at 1 kpc from the plane.
In addition to the large-scale variations, the median azimuthal velocity shows small-amplitude (a few km s −1 ) variations with galactic radius, with maxima at R ∼ 6.5 kpc (for φ > 0 and Z in [−600, 200] pc), R ∼ 8 kpc (for φ > 0 and Z in [−600, 600] pc), and R ∼ 10 kpc (Z in [−600, −200] pc). In the face-on maps (Fig. C.3), in which the colour range dynamics was reduced to heighten the contrast between velocity features, these maxima are visible as red circular arcs. Super-imposed on this large-scale variation, the azimuthal velocity also shows arc-shaped oscillations with small amplitude on a kiloparsec scale. Figure 14 shows a global increase in median vertical velocity, from the inner to the outer disc, but with complex vertical and azimuthal dependencies. The face-on ( Fig. 10 lower left panel and Fig. C.5) and edge-on maps (Fig. 11 lower left panel and Fig. C.7) show kiloparsec large, negative (green to blue) and positive (light green to red) velocity features, with an elaborate 3D geometry. Figure 15 presents the vertical projection ofṼ Z as a function of height Z for different azimuth slices and ranges in Galactic radius. In the outer disc (R > 10 kpc), the positive velocity feature appears inclined with respect to the Galactic plane, that is, it is located below the mid-plane at φ −15 deg, extending over most of the width of the plane for φ in ∼[−15, +15] deg and located mainly above the plane for φ 15 degrees. Still in the outer disc and for φ ∈ [−15, +15] degrees, the median vertical velocity is mildly symmetric with respect to the mid-plane, with a minimum at around  The bending velocity is negative (i.e. oriented towards the south Galactic pole) at negative azimuth for |Z| ∈ [0, 400] pc and in the inner disc at larger distance from the mid-plane. It is positive in the outer disc. Close to the Galactic mid-plane, the signal is weak and localised. It becomes stronger and spatially more extended with greater distance from the mid-plane. The absolute value of the breathing velocity is mostly lower than 1 km s −1 for |Z| < 800 pc. In the range |Z| ∈ [800, 1200] pc, the breathing velocity is partly positive in the first, second, and fourth quadrants, and it is negative in the third.

Vertical velocity
Using SEGUE spectra, Widrow et al. (2012) studied the vertical variations in mean vertical velocity,V Z , of a sample of high Galactic latitude (|b| ∈ [54, 68] degrees) outer disc stars (Galactic longitude l ∈ [100, 160] degrees). The mean vertical velocities they measured show a vertical asymmetry, with V Z < 0 km s −1 below ∼0.5 kpc and positive above. The mean vertical velocity also presents some oscillations. In the following year, Williams et al. (2013) studied the velocity field in an area of about 2 kpc around the Sun. Their (R, Z) maps show inversions of the sense of the mean vertical motion of the stars along the Z axis, producing zones of compression and zones of rarefactions of the stars. Recently, Carrillo et al. (2018) compared the velocity field derived with different proper motion catalogues and found great differences in particular in the median vertical velocity,Ṽ Z , maps. With the Gaia DR1 TGAS catalogue, they observed a breathing mode (a median motion of the stars away from the plane) in the inner disc and a downward bending beyond the Sun, over a distance of about 1 kpc.
The complex radial, vertical, and azimuthal dependencies of the vertical velocity make a comparison of samples with different selection functions difficult. The stars selected at positive azimuth and less than 2 kiloparsecs beyond the Sun (orange curves in the lower panels of Fig. 15) have some intersect with the sample of Widrow et al. (2012). Although not identical, the vertical velocity profiles look compatible. In the inner disc and φ > −15 degrees, we observe an increase in vertical velocity, with Z having similarities with the vertical profile of the RAVE-TGAS sample 7 of Carrillo et al. (2018), but with smaller amplitudes at large Z and a less pronounced symmetry with respect to the mid-plane (our inner discṼ Z are mostly negative). It should be noted that because the median V Z values are relatively modest, a small change in the vertical velocity zeropoint can modify the position of the inversion of the vertical motion. Figures 16-18 show the dispersions of the three galactocentric components of the velocities, σ V R , σ V φ , and σ V Z , as a function of galactic radius for negative (left) and positive (right) azimuths and for different Z layers (the different curves). The three velocity dispersions decrease with increasing radius. The gradient is significantly stronger at intermediate and large Z than in the midplane, with the vertical velocity dispersion σ V Z showing almost no gradient in the Z layer [−200, 200] pc. The dispersions are very symmetric with respect to the Galactic mid-plane, with the curves of symmetric layers showing very similar behaviours, including some kiloparsec-scale bumps/oscillations.

Radial, azimuthal, and vertical velocity dispersions
As shown on the right side of Fig. 11, the iso-velocity dispersions flare outwards. Two effects can act together to produce these flares. On the one hand, there is a radial evolution in the relative proportion of the short-scale length thick disc and the colder longer-scale length thin disc. On the other hand, 7 and distances from Astraatmadja & Bailer-Jones (2016).  with increasing outward distance, the vertical component of the gravitational force weakens, and for the same velocity, a star can reach larger distances from the mid-plane.
The velocity dispersions, in particular σ V R and σ V φ , show small-amplitude fluctuations that extend on a kiloparsec scale both radially and vertically. The face-on view of the disc (Fig. 10) shows that these hot features have a semi-circular geometry that extends at least 20 to 30 degrees in azimuth.

Discussion
The Milky Way is not an axisymmetric system at equilibrium. In the past few years (less than a decade), asymmetric motions (Casetti-Dinescu et al. 2011), gradients (Siebert et al. 2011), and wave patterns (Widrow et al. 2012) have been detected in the velocity field and were studied in increasingly more detail (Williams et al. 2013;Carlin et al. 2013;Sun et al. 2015;Carrillo et al. 2018;Pearl et al. 2017;Tian et al. 2017;Liu et al. 2017;Baba et al. 2018;Schönrich & Dehnen 2018). The second Gaia data release now offers a full-sky 3Dview of the complex Milky Way velocity pattern. It shows streaming motions in all three velocity components as well as small-amplitude fluctuations in the velocity dispersions.
Streaming motions might be produced by internal mechanisms (e.g. response of the stars to the bar and/or spiral structure) or by external perturbers (e.g. satellite accretion(s), impact of low-mass dark matter halos), or by combinations of both. It is beyond the scope of this paper to model the observations in detail. Below, we review some of the results of the already rich literature and discuss them with regard to the Gaia DR2 maps. Siebert et al. (2012) compared different models with two and four long-lived spiral arms. They successfully reproduced the gradient they had found the year before (Siebert et al. 2011) with a two-arm model and a specific location of the Sun, near the inner ultra-harmonic 4:1 resonance. Monari et al. (2014) showed that the bar can also induce a negative radial velocity gradient. Faure et al. (2014) studied the stellar velocity response to stable spiral perturbations. Within corotation (located in their model at 12 kpc), the model produces inward motions within the arms and outward motions between the arms. The model also induces vertical breathing modes, with stars moving away from the plane at the outer edges of the arms and towards the plane at the inner edges (still within corotation). Debattista (2014) also obtained breathing modes, with compression where the stars enter the spiral arms, and expansion where they exit. Monari et al. (2016a) developed an analytical model, based on phase-space distribution functions, to study the perturbations induced by a spiral potential. The model predicts breathing modes. Grand et al. (2016) used cosmological simulations to study the large-scale motions induced by the spiral arms in a Milky Way-like galaxy. The simulation shows radially outwards and azimuthally backwards motions on the trailing edge of the arms, while on the leading edge, the effect is reversed: the streaming motion is oriented inwards and forwards (see also Antoja et al. 2016). Monari et al. (2016b) studied the combined influence of the bar and two quasi-static spiral arms. The model produces horizontal (i.e. radial and azimuthal) streaming motion dominated by the influence of the bar and vertical breathing modes with spiral arms shape, but with the bar heightening the amplitude of the modes and shifting their locations. The vertical waves produced by internal mechanism models are usually breathing modes. Using N-body simulations, however, Chequers & Widrow (2017) recently showed that even in isolated Milky Way-like galaxies, random noise in the distributions of halo and bulge stars can produce long-lived bending waves in the disc that are observable beyond the solar circle. Figures 19 and 20 show the face-on maps of the median radial and azimuthal velocities, respectively, for the mid-plane layer ([−200, +200] pc). Two models of spiral arms are overplotted. The two-arm model of Drimmel (2000), derived from near infra-red data, is represented by thick black lines, and the locus of the minimum density between the two arms is shown by the thick dashed line. The spiral arms model of Reid et al. (2014) is represented with thin colour-coded lines (see caption of Fig. 19). Reid et al. (2014) used masers as tracers of the spiral arms. It should be noted that these masers are associated with massive stars that are much younger than the giant stars whose kinematics is mapped in this section. The Local Arm shows some coincidence with the ridge of negative median radial velocities, and its trailing edge is close to the boundary between positive and negativeṼ R . This might even be fortuitous as the Local Arm is usually considered a weak structure (Churchwell et al. 2009). The locus of minimum density between the two near-infrared arms also matches the boundary between positive and negative median radial velocities. The locus also correspond mildly with the semi-circular faster azimuthal velocity pattern (yellow-red arc in Fig. 20). Dynamical models of bar and/or spiral arms predict streaming motions and changes in sign of the median radial velocity. It is therefore tempting to see a link between the radial velocity oscillation and the near-infrared arms. Siebert et al. (2012) indeed reproduced the negative radial gradient with a two-arm model. On the other hand, it should also Fig. 19. Face-on map of the median radial velocity (in km s −1 ) for the mid-plane layer ([−200, +200] pc), derived using the giant sample. The two-arm model of Drimmel (2000), adjusted on near infra-red data, is over-plotted as thick black lines. The thick dashed line highlights the locus of minimum density between the two arms. The spiral arms model of Reid et al. (2014) is also over-plotted, i.e. from the inner to the outer disc: Scutum (cyan), Sagittarius (magenta), Local Arm (blue), Perseus (black), and the Outer Arm (red).  Fig. 19 for the median azimuthal velocity,Ṽ φ . be noted that Liu et al. (2017) obtained a radial oscillation by adjusting the positive radial gradient with a bar model. The mapping in 3D ofṼ R andṼ φ brings new constraints to the models.
Vertical to the Galactic disc, we expect the kinematics to reflect the large-scale warp. If the Milky Way warp is a long-lived structure, then we expect an associated kinematic signature towards the Galactic anti-centre in the vertical velocities. Figure 10 (lower left plot) indeed seems to exhibit a systematic vertical velocity of about 2-3 km s −1 at R = 10−11 kpc in the direction of the anti-centre. However, this signal is weaker than expected from current empirical descriptions of the stellar warp, which assume the warp to be stable and non-precessing, and might indicate that the warp is instead an unstable transient feature. It is worth comparing theṼ Z map of the giants (Fig. 10, lower left panel) to an equivalent map for the young OB sample. Lineof-sight velocities are not available for this sample (see Sect. 2.1 and Katz et al. 2018 for details), so that we cannot directly calculate V Z for each star. However, at low Galactic latitudes, we can estimate a V Z since most of the vertical motion is seen in the proper motions perpendicular to the Galactic plane. The vertical velocity where S = U cos l + V sin l, and similarly for S , which contains both differential Galactic rotation and the own peculiar motion of the star parallel to the Galactic plane, S * . Neglecting the latter and assuming a flat rotation curve, we estimate V Z by taking S ≈ V LS R (R /R−1) sin l in the above equation. Using stars in the OB sample within 200 pc of the Galactic plane, we map the median V Z , shown in Fig. 21 (our approximation means that we have effectively ignored a S * tan b for each cell). The resulting map is distinctly different from the map for the giants, and show no signs of a warp signature. We recall that our sample traces the motions of the gas from which these stars have recently been born. The lack of any warp signature here therefore again suggests that the warp is an unstable transient feature.
External perturbations by satellites or dark-matter haloes, for instance, can also disturb the disc structure and velocity field, and may in particular induce vertical waves (Widrow et al. 2014;Feldmann & Spolyar 2015). Simulating the interaction between a satellite and a galactic disc, Widrow et al. (2014) showed that the response of the disc depends, in particular, on the relative vertical velocity of the satellite. A satellite approaching "slowly" will primarily bend the disc, while a fast satellite will produce breathing and higher order modes. The Sagittarius dwarf galaxy has been interacting with the Milky Way for several Gyr. It is therefore a natural suspect for the perturbations observed in the velocities as well as in the outer disc. Simulations of the accretion of Sagittarius by the Milky Way were indeed able to reproduce the local vertical velocity pattern measured by Widrow et al. (2012) Laporte et al. 2018b), but also large-scale outer-disc features such as the Monoceros ring (Laporte et al. 2018b). The Large Magellanic Cloud can also induce vertical modes in the local disc, but with significantly smaller amplitudes than Sagittarius (Laporte et al. 2018b,a). External perturbers do not only modify the vertical velocities, but also the horizontal ones. The effect has in particular been studied in the velocity plane Gómez et al. 2012b) and is discussed in the next section.
The vertical velocity field (Sect. 3.4) is quite complex, with different behaviour in the inner and outer disc and radial, azimuthal, and vertical dependencies. As previously observed by Carrillo et al. (2018), it cannot be described by a single bending, breathing, or higher mode. We likely witness a superposition of modes that may be of several different origins.

Kinematic substructure
In this section, we revisit the substructure in the velocity plane of the solar neighbourhood and explore it also in the velocity distribution of distant regions from the Sun. We focus on the in-plane velocities because they have been demonstrated to show most of the substructure. We remark that a detailed study of the kinematic substructure present in local and distant neighbourhoods is beyond the scope of the present study. Here we present only a first exploratory look and focus on the quality of the Gaia data, the aspects that allow us to scientifically verify the data, and the highlights of our findings.

Kinematic substructure in the solar neighbourhood
As explained in Sect. 2.4, we built our solar neighbourhood sample by selecting stars closer than 200 pc to the Sun. The number of stars (366 182) is an order of magnitude larger than local samples such as Geneva-Copenhagen Survey (GCS, Nordström et al. 2004) or RAVE (Steinmetz et al. 2006). Moreover, the uncertainties in velocity of this sample, which are smaller than 1 km s −1 for about 80% of the sample, allow us to probe substructure scales that are significantly smaller than ever before. Figure 22 shows a 2D histogram with a bin size of 1 km s −1 of the in-plane Cartesian heliocentric velocities U and V of the solar neighbourhood sample. This velocity distribution is highly structured. We observe many nearly horizontal arch-like structures that have never been seen before. Even the dynamical stream of Hercules, located at negative U velocities and V −50 km s −1 , now appears to be split into at least two of these branches (at V −38 km s −1 and V −50 km s −1 ) and perhaps the stream at V −70 km s −1 is also associated with the same structure. These arches appear for the whole range of V. We note, for instance, the new low-density arch at V 40 km s −1 . Some of them are not centred on U and others are inclined in V. Additionally, there is a clear under-density of stars also with an arched shape that extends from (U, V) (−100, −25) to (U, V) (75, −65) km s −1 immediately above the Hercules stream, which separates the velocity plane in two. This gap is not as horizontal as the over-dense arches. Figure 22 also shows more strongly rounded structures with sizes of about 10 km s −1 , especially at the centre of the distribution, which correspond to previously known moving groups and dynamical streams. We also see small high-density clumps that might be associated with known open clusters in the vicinity of the Sun. For instance, the yellow clumps at (U, V) = (−42, −19) and (U, V) = (−7, −28) km s −1 correspond to the Hyades and Pleiades clusters, respectively. All of these substructures, that is, the medium and small structures, appear to be embedded in the larger arched substructures.
To study the substructure at different scales, we have performed the wavelet transform of the velocity plane of our neighbourhood sample. This is a mathematical transform that decomposes an image into the basis of the so-called wavelet function of different sizes. In practice, this transform yields a series of planes containing the substructures at different scales, and it has been extensively used to analyse the velocity substructure (e.g. Antoja et al. 2008;Chereul et al. 1999). A thorough description of this technique can be found in Starck & Murtagh (2002).
In Fig.23 (top) we show the wavelet transform at scales of 8 km s −1 . This figure looks very similar to the wavelet transform applied to previous local sample (e.g. Antoja et al. 2011), with several rounded structures organised in a larger arched structure. The green and orange crosses mark peaks that are >3σ and 2−3σ significant, respectively. Most of the peaks detected at this scale are already known moving groups and dynamical streams. For comparison, we superpose the peaks detected with RAVE in Antoja et al. (2011) with squares following the same colour code. The most prominent structures are the moving groups of the Hyades ((U, V) = (−33, −16) km s −1 ) and Pleiades ((U, V) = (−11, −24) km s −1 ), which must not be confused with the open clusters with the same name, and the Sirius moving group (U, V) = (10, 3) km s −1 . There are some groups that now have a larger significance than before, such as the two peaks inside the known Hercules stream (two green crosses at V −50 km s −1 ). Moreover, new significant groups appear, such as the substructures at (U, V) (−90, 0) km s −1 .
The kinematic substructure of smaller size can also be explored thanks to the precision of our data. We show it here in Fig. 23 (bottom) with the wavelet transform at scales of 2 km s −1 . The plot shows a zoom-in into the central parts of the velocity plane. The figure again reveals prominent arched over-densities of stars, each one with a different inclination. We distinguish here at least five of these arches. Each of them shows internal smaller and rounder substructures, some of which mark the peaks of nearby clusters.
The red circles in this plot mark the velocity positions of eight nearby well-defined and fairly rich open clusters that are located at distances closer than 200 pc (Hyades, ComaBer, Pleiades, Praesepe, alphaPer, IC2391, IC2602, and NGC2451A). Their heliocentric velocities have been computed from mean parameters provided in Babusiaux et al. (2018). The Hyades and Pleiades clusters have the largest number of stars in the sample, thus showing high-intensity peaks in the wavelet transform at these small scales. Other clusters do no show a particular overdensity in the wavelet space because they might be composed of a smaller number of stars in the local sample. Almost all clusters fall inside the larger substructure that is formed by the arches, and seven out of eight lie in the same arch.

Kinematic substructure in distant regions
To study the distant velocity distributions, we have selected stars from the main sample that are located in the Galactic disc with a cut of |Z| < 400 pc and focused on a ring around the Sun between distances on the Galactic plane of 0.5 and 1 kpc. We furthermore partitioned this ring into sectors of 45 deg centred at the longitude of l = [0, 45,90,135,180,225,270,315] deg. The number of stars in these regions ranges from 90 831 to 171. We then used Galactic cylindrical velocities since they are in a more natural reference system for these distant regions. We plot −V R instead of V R to obtain the same orientation as U (i.e. velocity towards the Galactic centre). The median The green and orange symbols mark peaks that are >3σ and 2 − 3σ significant, respectively, for DR2 (crosses) or for the RAVE sample (squares) from the study of (Antoja et al. 2012; their table 2). Bottom panel: scales of around 2 km s −1 . The green crosses mark peaks that are >3σ significant. The red circles mark peaks for the velocity positions of nearby open clusters that are located at distances closer than 200 pc (see text). This panel shows a zoom-in on the central parts of the velocity plane. velocity uncertainties in the whole ring are of ( V R , V φ , V Z ) = (1.2, 1.3, 0.6) km s −1 , and 80% of the stars have uncertainties in all velocity components smaller than 3. km s −1 .
The 2D histograms of these eight distant regions are shown in Fig. 24. In contrast to Fig. 22, we now use a bin of 2 km s −1 to account for the larger uncertainty in these samples. The panels are oriented such that the Galactic centre would be located at the left-hand side of the figure. For comparison, we include the distribution of the solar neighbourhood in the middle panel with the same bin size.
Substructure is ubiquitous in these panels. Although we note a loss of definition in the substructures compared to the local volume because of the larger uncertainties, we still see clearly arched structure. Additionally, we note that the velocity structures change from region to region, with greater changes in Galactic radius than in azimuth (i.e. greater changes for the different columns than for the different rows). For example, the gap that separates the Hercules stream clearly moves from smaller to larger V φ from the outer regions (three rightmost panels) through the regions at the solar Galactic radius (three middle panels) to the inner regions (three leftmost panels). Effectively, the gap moves from V φ 200 km s −1 to V φ 240 km s −1 over a distance range of 2 kpc. Although these findings are consistent with previous studies Monari et al. 2017), the resolution is unprecedented and we did not have to treat the data with any sophisticated method, as was required in previous work.
In most of the panels, but especially at the rightmost ones, this resolution allows us for the first time to see a structure below the Hercules stream that is separated by a secondary gap. This gap is located at around V φ 200 km s −1 in the three panels on the right, but lies at a different velocity in the other panels. In addition, the uppermost arch (at the highest V φ ) that is observed in the solar volume is located at even larger V φ in the inner Galactic regions (three left panels) and at smaller V φ in the outer quadrants.

Discussion
We report a new global arrangement of stars in the velocity plane of the solar neighbourhood in which stars are organised in arches with nearly constant heliocentric velocity V that extends to a wide range in U. This discovery is made possible by the higher accuracy and greater number of stars in the Gaia local sample, which exceeds the accuracy and number of stars of previous data sets.
Similar kinematic arches have appeared in several simulations, for example in the models by Dehnen (2000) and Fux (2001), who both studied the effects of the Galactic bar on the local velocity distribution. The velocity arches appeared in their simulations, which use backwards integration, and the number of arches depended on the total integration time. Fux (2001) related these arches to phase-mixing. Similar to what we observe here, their arches are not perfectly centred on U or inclined in V. Similarly, Antoja et al. (2009) explained that arches also appeared in forward integrations in axisymmetric potentials as a result of a non-relaxed population and that these arches were not centred on U when the bar was added to the potential.
Arched kinematic substructures have also been reported in simulations of spiral arm effects in Quillen & Minchev (2005) and Antoja et al. (2011), and structures elongated in U were observed in the simulations of De Simone et al. (2004), with transient episodes of spiral structure, although in the latter case they did not appear as strongly curved as in the Gaia data. The scenario of transient spiral arms with the superposition of multiple waves with different frequencies is independently supported by several other pieces of evidence, especially from simulations (see e.g. Sellwood & Carlberg 2014). Thus, the possibility that spiral arms or the possible transient multiple waves generate one or several of these arches should be explored in future work as well, in addition to the possibility that effects of phase-mixing might also play a role because these modes are transient.
Clear arches are also observed in the models by Minchev et al. (2009) and have been followed by Gómez et al. (2012b) in a phenomenon that they dubbed "ringing", which these authors attribute to the impact of a satellite on the Galaxy disc. Using orbital integrations, a semi-analytical method, and N-body simulations, these authors found that as a response to the perturbation, the disc experiences a series of waves of constant energy that propagate through the disc and manifest themselves as arches in the U − V plane. All these authors explained that these arches can give information on the time of the perturbation, the orbital parameters of the perturber, and its mass. It seems therefore that the arches are evidence of ongoing phase-mixing Fig. 24. Velocity plane of the stars in several distant regions. The eight outer panels correspond to heliocentric sectors in the Milky Way disc selected from a ring between 0.5 and 1 kpc from the Sun, in a layer of |Z| < 400 pc and 45 deg of amplitude centred at the Galactic longitudes indicated in each panel. The panels are 2D histograms with a bin size of 2 km s −1 , thus, the colour scale indicates the number of stars per 4 ( km s −1 ) 2 . The middle panel shows the solar neighbourhood sample (defined in Sect. 2.4) for comparison (the same as in Fig. 22, but with a larger bin size). The number of stars is indicated in each panel. We have reversed the horizontal axis V R to obtain the same orientation as U (i.e. velocity towards the Galactic centre).
in the disc, but this might be attributed to several causes such as a perturbing satellite, the formation of the bar, or other changes in the Galactic potential to which the stellar disc is still adapting. Minchev et al. (2009) searched for observational signatures of this scenario and focused on thick-disc samples, arguing that these stars do not experience additional perturbations by molecular clouds or spiral arms and that these arches might therefore be better distinguished. These authors suggested that some of the previously known streams might be signatures of a minor merger. However, because of the low number of stars and the velocity resolution at that time, these structures did not have a clear arched shape, and as the authors admit, were not significant enough. Recently, Gómez et al. (2012a) found a few significant peaks in the energy distribution of the Schuster et al. (2006) and Lee et al. (2011) catalogues that might well be linked to ringing and that showed rough similarities with those produced by a model of the Sagittarius dwarf galaxy.
For the first time, the arched substructures are clearly observed here, and surprisingly, they have appeared in a sample that is dominated by thin disc stars, since no selection on metallicity has been made here. We note their different inclinations and ranges in U, which must contain information on the mechanism that causes them. In any case, a detailed comparison with the different models is required to ascertain which exactly are the more plausible perturbing agents that may have caused a non-equilibrium state in the disc.
At larger scales, the velocity distribution that Gaia has measured shows rounded structures that are perfectly consistent with previous studies, and that might be of dynamical origin, related to the spiral arms and bar gravitational potential, as suggested originally by Kalnajs (1991) and later by many others (see Sect. 1). It is also worth noting that the open clusters considered here, which are more prominent inside the solar neighbourhood, are all located on top of the arches in the velocity plane and thus seem to participate in a common dynamics. We note, however, that this is not a complete sample of clusters and that further study is required to extract definitive conclusions. If it is confirmed that most of the clusters follow this arched organisation, the ages of the clusters might give us information on the timescales of the agent that created these arches. The clusters studied here are quite young, which might indicate that the causing agent(s) have been acting until quite recently.
The Gaia data up to distances of 1 kpc in the Galactic plane show that the arches also exist in more distant regions. The quality of the data will certainly allow us to quantify the changes in substructure with position. We already observe movement in the velocity of the gap between the Hercules stream and the remaining distribution as a function of radius, specifically at lower V φ for larger Galactic radius. This was first discovered by Antoja et al. (2014) with RAVE data, based on which the authors demonstrated that this change is compatible with the orbital effects of the Galactic bar near the outer Lindblad resonance. This finding was also confirmed later on by Monari et al. (2017) using TGAS and LAMOST. However, while in Antoja et al. (2014) a special projection of the velocities and complex analysis tools were required, here the gap and its variation with position is unambiguous by a simple inspection of data. We also note that the range in Galactic radius that can be explored is larger in the present study (2 kpc) than in these previous studies ( 1 kpc). Hunt et al. (2018) also showed some variations in Hercules as a function of distance with APOGEE-2 South data, but only for a specific range of Galactic radius, because their data were limited to a single line of sight. The thorough study of this velocity variation with Gaia data will lead to a better determination of the bar properties, such as pattern speed and orientation angle.

Conclusions
Taking advantage of the full-sky coverage, precise distances, and velocities of Gaia DR2 for more than three million giant stars, we have mapped the velocity field of the galaxy in 3D over a large portion of the disc (5 ≤ R ≤ 13 kpc, −30 • ≤ φ ≤ +30 • , and −2 ≤ Z ≤ +2 kpc). The picture of the Milky Way disc kinematics drawn by Gaia DR2 is both rich and complex. Streaming motions are observed in all three velocity components, andṼ Z likely shows a superposition of modes. The velocity dispersions also show some small-amplitude perturbations superimposed on a large-scale radial decrease that is quite symmetric in azimuth and with respect to the Galactic mid-plane. This is also the first time that the velocity precision is high enough and the statistics is large enough to resolve the small scales of the velocity plane of the solar neighbourhood. We find that stars clearly appear to be organised in kinematic arches that are oriented approximately along the horizontal direction in the U − V plane. The origin of these arches is probably related to the non-equilibrium of the Galactic disc, which could have been induced by several causes or by a mix of them. For instance, similar arches have been seen in simulations of a satellite impact on the disc. These arches are also observed in the Gaia DR2 data of more distant regions, and we also see the change in some of the previously known kinematic groups as a function of position in the disc in a manner consistent with models of the orbital effects of the Galactic bar. The velocity plane of the local neighbourhood and surroundings certainly appears to be as complex as ever, with a variety of structures of many different scales and shapes. Beyond these substructures, which are signatures of current and past events of the formation and evolution of the Milky Way, this first inspection of the Gaia DR2 data clearly promises exciting discoveries and significant advances in this field in the coming years.
Gaia DR2 is now available to the astronomical community so that a variety of questions about the Milky Way can be addressed that far exceed the specific examples visited in this paper. In particular, Gaia DR2 provides transverse velocities for several hundred million stars. Although it was beyond the scope of this paper to use them, there is no doubt that they have much tell about the Galaxy kinematics and dynamics. Moreover, the joint use of Gaia DR2 and large ground-based spectroscopic surveys will allow combining the precise Gaia DR2 distances and velocities with detailed chemistry. We can expect this great wealth of information to trigger an intense activity in the galactic community in the years to come. Following Gaia DR2, Gaia DR3 is scheduled to be released in a few years 8 . The many promises of the next release include a further improvement in the precision and accuracy of the astrometric, photometric, and spectroscopic data, an increase by a factor 5 to 10 of the stars with full velocities, detection and characterisation of multiple systems, and chemistry for millions of stars.
The lower and upper 1σ uncertainties on the estimations of the median velocities are calculated as using the same notation as above.
The lower and upper 1σ uncertainties on the estimation of the velocity dispersions are calculated as   for an additional selection of stars with velocity errors in each component smaller than 2 km s −1 , respectively. It is evident from the comparison of all panels in the two figures that this cut in velocity error introduces strong biases in the median velocity field. The differences in median velocity can be of up to 20 km s −1 , and the global appearance of the velocity field changes substantially. The biases result from the correlations between velocities and velocity uncertainties. Selecting on the velocity uncertainties modifies the shape of the velocity distributions and therefore biases the measures of the moments of these distributions.