3D kinematics and age distribution of the Open Cluster population

Open Clusters (OCs) can trace with a great accuracy the evolution of the Galactic disk. The aim of this work is to study the kinematical behavior of the OC population over time. We take advantage of the latest age determinations of OCs to investigate the correlations of the 6D phase space coordinates and orbital properties with age. We also investigate the rotation curve of the Milky Way traced by OCs and we compare it to that of other observational or theoretical studies. We gathered nearly 30000 Radial Velocity (RV) measurements of OC members from both Gaia-RVS data and ground based surveys and catalogues. We computed the weighted mean RV, Galactic velocities and orbital parameters of 1382 OCs. We investigated their distributions as a function of age, and by comparison to field stars. We provide the largest RV catalogue available for OCs, half of it based on at least 3 members. Compared to field stars, we note that OCs are not exactly on the same arches in the radial-azimuthal velocity plane, while they seem to follow the same diagonal ridges in the Galactic radial distribution of azimuthal velocities. Velocity ellipsoids in different age bins all show a clear anisotropy. The heating rate of the OC population is similar to that of field stars for the radial and azimuthal components but significantly lower for the vertical component. The rotation curve drawn by our sample of clusters shows several dips, which match the wiggles derived from non-axisymmetric models of the Galaxy. From the computation of orbits, we obtain a clear dependence of the maximum height and eccentricity with age. Finally, the orbital characteristics of the sample of clusters as shown by the action variables, follow the distribution of field stars. The additional age information of the clusters points towards some (weak) age dependence of the known moving groups.


Introduction
The study of kinematical properties of the Milky Way open clusters (OCs) has a long tradition. Their motion can be used to understand the Milky Way's gravitational potential and the various perturbations which act on the structure and dynamics of the Galactic disc. The solar neighbourhood is known to be clumpy (e.g. Eggen 1996;Dehnen & Binney 1998;Antoja et al. 2012), and the relation of observed substructures with some OCs, in particular the Hyades, has been established (Eggen 1958;Chereul et al. 1999). While stellar streams were initially believed to be remnants of star clusters, their origin is now believed to be the result of the disruption of star clusters or of dynamical origin such as the resonant trapping by the bar The tables with star and cluster velocities are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/?/? and the spiral arms (Famaey et al. 2008) or the passage of the Sagittarius dwarf galaxy within the galactic plane (Monari et al. 2018;Khanna et al. 2019). Stellar streams could also originate from a combination of the aforementioned processes that happen together. Gaia DR2 data (Gaia Collaboration et al. 2018b) revealed the complexity of the local velocity distribution with an unprecedented resolution, showing in particular that the velocity distribution of nearby OCs overlaps well with prominent arched over-densities of moving groups (Gaia Collaboration et al. 2018d).
Gaia DR2 has considerably fostered the study of OCs with the determination of new memberships for an unprecedented number of stars and clusters. A first large catalogue was published by Cantat-Gaudin et al. (2018) who systematically looked for members around the ∼3300 catalogued OCs (mostly from Dias et al. (2002) and Kharchenko et al. (2013) studies), only Article number, page 1 of 16 arXiv:2012.04017v1 [astro-ph.GA] 7 Dec 2020 A&A proofs: manuscript no. main based on Gaia DR2 positions, parallaxes and proper motions. Membership probabilities were computed for ∼400 000 stars and among other parameters, the most probable distances were determined for 1229 OCs. A large fraction of candidate OCs could not be confirmed by Gaia or was proved to correspond to chance alignment. The number of OC candidates has also significantly increased, either thanks to serendipitous discoveries in Gaia DR2 (Cantat-Gaudin et al. 2018;Ferreira et al. 2019), or as the result of systematic searches (Castro-Ginard et al. 2018, 2019Cantat-Gaudin et al. 2019;Sim et al. 2019;Liu & Pang 2019). An updated catalogue of membership probabilities for 1481 OCs has been provided by . The physical and kinematical properties of OCs has also been revisited thanks to Gaia. Using the memberships from Cantat-Gaudin et al. (2018) and only Gaia data, Bossini et al. (2019) determined the age, distance modulus and extinction of 269 OCs, and Soubiran et al. (2018a) computed the 6D phasespace information of 861 clusters. Carrera et al. (2019) increased by 145 the number of OCs with full 6D phase-space information by looking for members in the GALAH and APOGEE spectroscopic surveys. The vertical distribution of young clusters was found to be very flat, with a dispersion of vertical velocities of 5 km s −1 , while clusters older than 1 Gyr span distances to the Galactic plane of up to 1 kpc with a vertical velocity dispersion of 14 km s −1 , typical of the thin disc. Recently,  derived physical properties of all known OCs identified in the Gaia data in a homogeneous fashion, with a method based on isochrone fitting and an artificial neural network, ending up with 1867 clusters with reliable ages.
The paper is organized as follows. Sect. 2 describes the collection of individual RVs gathered from Gaia DR2 (Sartoretti et al. 2018;Katz et al. 2019) and from ground-based resources, used to build two catalogues of mean RVs, per star and per cluster. In Sect. 3, 6D galactocentric coordinates are computed. We compare the resulting distributions with those of nearby field stars. We evaluate how the different coordinates vary with the age of the clusters and we determine the age velocity relation of the population. We investigate how OCs follow the theoretical rotation curve of the Galactic disc. Orbital parameters and actions are described in Sect. 4 and used to revisit the link between OCs of different ages and the phase-space substructures of the disc. The conclusions of this study are summarized in Sect. 5

Input data
We took advantage of the catalogue of OCs by  which provides a list of probable members for 2017 OCs that they used to estimate ages. Their memberships mostly come from Cantat-Gaudin & Anders (2020) who used the unsupervised classification scheme UPMASK (Krone-Martins & Moitinho 2014;Cantat-Gaudin et al. 2018). They also applied UPMASK to the recently discovered University Barcelona Clusters (UBC, Castro-Ginard et al. 2018, 2019 and to those discovered by Liu & Pang (2019). For the Hyades and Coma Berenices, they adopted the list of members published by Gaia Collaboration et al. (2018a), UPMASK being unable to recover members for populated clusters that are too extended on the sky.
Based on the aforementioned list of members, containing ∼475 000 stars, we gathered all the RV measurements available for them in various surveys and catalogues. Our main source was Gaia DR2 which includes RVs for about 7 million stars (Sartoretti et al. 2018;Katz et al. 2019). We also queried sev-eral catalogues from large spectroscopic surveys in addition to Gaia: the latest public version of Gaia-ESO survey (Randich et al. 2013), APOGEE DR16 (Ahumada et al. 2020), RAVE DR6 (Steinmetz et al. 2020), and GALAH DR3 (Buder et al. 2020). We did not consider LAMOST (Cui et al. 2012;Xiang et al. 2015) because its RV precision and accuracy (∼5 km s −1 ) are not at the same level. We included radial velocities derived in by the OCCASO survey (Casamiquela et al. 2016, Carrera et al. in preparation). We also considered the RV catalogues by Soubiran et al. (2018b); Mermilliod et al. (2009Mermilliod et al. ( , 2008; Worley et al. (2012); Nordström et al. (2004). Several quality cuts were applied to individual measurements. For RAVE we applied the criteria suggested by Steinmetz et al. (2020) by selecting stars verifying |correctionRV| <10 km s −1 , σ RV < 8 km s −1 and correlationCoeff >10. For APOGEE, we rejected the stars flagged with VERY_BRIGHT_NEIGHBOR, BAD_PIXELS or LOW_SNR as recommended in the online documentation of the survey. For Gaia RVS we filtered the erroneous RVs found by Boubert et al. (2019). Despite these cuts, there were still individual measurements with large uncertainties incompatible with the required precision for this work. This convinced us to filter out all the individual RVs having uncertainties larger than 8 km s −1 , the same cut as for RAVE. The rejected values represent 6% of the full set, 10% for the RVS set. Then we noticed 43 stars with |RV| > 200 km s −1 which turned out to be mainly OB-type or Wolf-Rayet or spectroscopic binary stars according to Simbad. Such stars with unreliable or variable RVs were rejected. Table 1gives the number of cluster members retrieved in each of the catalogues after filtering, with the median uncertainty of the corresponding RVs, as quoted in the catalogues. In Soubiran et al. (2018b), the uncertainty corresponds to the standard error of the weighted mean for stable stars that have been followed-up for exoplanet detection. The high precision of individual measurements and number of observations explain the very low median uncertainty of that catalogue (0.002 km s −1 ). APOGEE, OCCASO and the catalogues from Mermilliod et al. (2009Mermilliod et al. ( , 2008 and Nordström et al. (2004) also contain stars that have been observed several times in order to identify binaries. In that case, the RV uncertainty corresponds to the quadratic sum of the single measurement error with the scatter of the measurements. RAVE also has multiple observations for a small fraction of stars and provides the individual values with their uncertainty corresponding to the error of a single measurement. Finally the Gaia-ESO survey includes RVs of the same star obtained with different setups that give different measurement errors (Jackson et al. 2015).
This star sample is dominated by dwarfs and has a median effective temperature Teff of 4874 K (based on Teff provided in Gaia DR2), with more than 13% of the stars being hotter than 6250 K. We note that uncertainties of stars hotter than 6250 K are significantly larger (median 1.75 km s −1 ) than that of cooler stars (median 0.53 km s −1 ). This mainly reflects the lower precision of the RVs of rotating stars, which are frequent among AFtype stars in young OCs. Fig. 1 shows the distribution of the RV uncertainties of each catalogue. Note that uncertainties are not determined the same way in the different catalogues. The 28 371 RV measurements of OCs members that we gathered correspond to 23 453 different stars.

Zero points of RV catalogues
This sample allowed us to assess the consistency of RVs in the different catalogues. The comparison of RVs for stars in common in two catalogues gives an idea of potential offsets due to zero-   Mermilliod et al. (2008Mermilliod et al. ( , 2009; RAV: RAVE DR6; OCC : OCCASO Casamiquela et al. (2016); S18: Soubiran et al. (2018b); NOR: Nordström et al. (2004); WOR : Worley et al. (2012) point differences, together with their typical precision. Zeropoint differences are a result of the different observing modes, instrumental characteristics and calibration procedures of each instrument. They have to be taken into account when combining RVs of different origins. A subset of 3 116 stars have measurements in two or more catalogues. The RV difference between catalogues that have more than 20 stars in common is presented in Fig. 2, and the corresponding comparison among catalogues in Table 2. The RVs have a general good agreement with offsets smaller than 0.5 km s −1 and between 0.5 and 1.4 km s −1 for comparisons that involve RAVE or GALAH. The dispersions (measured by the median absolute deviation MAD) are typically of the order 1 km s −1 or below, consistent with the precision of the catalogues listed in Table 1 for all combinations of surveys. In these comparisons of the different catalogues, we did not see any significant trend with colour or apparent magnitude. Most of the RVs available for the cluster members are provided by Gaia DR2 or the Gaia ESO survey. GES has indeed a large fraction of its observing program dedicated to OCs and nicely increases the number of clusters for which a RV can be  Table 2. computed. Gaia and GES have an offset of ∼0.3 km s −1 , similar to the offset between GES and APOGEE, while Gaia and APOGEE agree at a level better than 0.1 km s −1 . The zero-point of RAVE is found to be different from the other surveys by ∼1 km s −1 (RAVE underestimates the RVs compared to the others) which is larger than what was reported in previous studies. Steinmetz et al. (2018Steinmetz et al. ( , 2020 and Sartoretti et al. (2018) report an offset of ∼0.3 km s −1 between RAVE and Gaia DR2. Here, even when we apply strong quality cuts on both RAVE and Gaia DR2, the median difference remains 1 km s −1 . RAVE also differs from APOGEE by 1.42 km s −1 , and from Mermilliod et al. (2008Mermilliod et al. ( , 2009 by 0.79 km s −1 . The agreement with S18 is, however, better at a level of 0.2 km s −1 , possibly related to the bright magnitude of the common stars and their colour range. We also note a systematic offset of GALAH DR3 with all the other catalogues, opposite to that of RAVE and of smaller amplitude, but still larger than the value of 0.22 km s −1 reported by Buder et al. (2020) for the comparison to Gaia DR2. These larger offsets are possibly related to the fact that our star sample is dominated by dwarfs, with a fraction of early-type stars larger than in RAVE and GALAH in general. For RAVE, we see indeed a trend of the offset being slightly larger for hot stars than for cool stars.
In order to put all the RVs on the scale of Gaia we applied a zero-point correction to the individual RVs from all the non-Gaia catalogues, according to the offsets listed in Table 2.

Mean radial velocities
In order to compute the mean RV of each OC, we first computed the mean RV of each star, since some of them (12%) have multiple measurements from different catalogues. For both the mean per star and the mean per cluster, we used a weighted procedure based on the errors of individual measurements, following Soubiran et al. (2013Soubiran et al. ( , 2018a. For each star or each cluster, the mean RV is computed by attributing to each RV measurement a weight ω i defined as ω i = 1/ 2 i , i being the RV error. Outliers were rejected on the basis of a 3σ clipping. A fraction of about 8% among stars with multiple measurements have a RV uncertainty larger than 3 km s −1 and we suspect that they are binaries with variable RV. Such stars are not rejected but the procedure gives them less weight when the mean RV of their parent cluster is computed. However, in the large majority of cases, we have only one measurement per star so that binaries cannot be identified. Binaries that have large variations of their RV may alter the mean RV of the parent OC in case of few members. The catalogue of stars provides the mean RV in the Gaia RVS scale, with its uncertainty and the number of measurements, as well as the membership probability from Cantat-Gaudin et al. (2020). It includes 23 424 unique stars, with 97 of them appearing twice owing to a non-null probability to belong to two different clusters being close to each other on the sky.
In order to compute the mean RV per OC, we considered only the stars with a membership probability higher than 0.4, the threshold value found by Soubiran et al. (2018a) to be the best compromise between the largest number of members and the lowest contamination by field stars. In the end, 1 382 OCs have a mean RV, which represents an improvement by 60% in size compared to the previous RV catalogue of OCs, based on Gaia only (Soubiran et al. 2018a). It also supersedes the catalogue by Carrera et al. (2019) based on Gaia, GALAH and APOGEE. Before Gaia DR2 the two largest compilations of cluster RVs were those of Kharchenko et al. (2013) and Dias et al. (2002) with respectively 962 and 703 objects. This is the first determination of mean RV for most of the recently discovered OCs. In particular half of the UBC clusters (Castro-Ginard et al. 2020) are part of our catalogue, and they represent nearly 20% of the full sample of OCs having a known RV. Among the 75 high confidence clusters recently discovered by Liu & Pang (2019), 35 are present in the latest catalog from Cantat-Gaudin et al. (2020) and we provide a RV for most of them. For instance, we provide the mean RV of UBC 274 (also reported as LP 5), a recently discovered disrupting OC (Castro-Ginard et al. 2020; Piatti 2020), based on 18 members observed with RAVE, GALAH and Gaia RVS (Fig. 3). Thanks to the combination of Gaia RVS and GES, several OCs have now their RV based on hundreds of members, the top two being Trumpler 5 and NGC 3532 with respectively 659 and 664 stars. For about 18% of the sample, the mean RV is based on more than 10 stars, and for 50% it is based on at least 3 stars. The median uncertainty of the weighted mean RV is 1.13 km s −1 when the full sample is considered.
Among the clusters with fewer members, some exhibit a large error which renders their mean RV uncertain. Some 430 OCs have a RV based on one star only, with no information on its potential variability (binarity). This represents 31% of our sample, nearly the same proportion as in Soubiran et al. (2018a). Selecting the most reliable OCs which have a RV uncertainty lower than 3 km s −1 based on at least 3 stars, we get 513 clusters. This sub-sample has a median uncertainty of 0.55 km s −1 and a me-  Table 2, and from Gaia RVS in orange, as a function of the G BP − G RP colour and G magnitude from Gaia. The grey line represents the weighted mean, RV=-23.05±0.25 km s −1 , based on 18 members. dian number of stars of 9. This is 107 more reliable OCs than in Soubiran et al. (2018a).
There are 21 high velocity OCs in the catalogue with |RV| > 100 km s −1 but 13 are based on a single star, therefore to be considered with caution. The 8 remaining high velocity clusters with more members were previously known.
The catalogue of OC velocities (released as online data) provides the RV per cluster with its uncertainty and the number of members used in the average, as well as galactic velocities, orbital parameters and actions computed in the next sections.

Open Clusters in phase space
In this section, we combined the mean OC positions, distances, proper motions from Cantat-Gaudin et al. (2020) and our mean RVs to compute heliocentric and Galactocentric cartesian and cylindrical positions and velocities for the full sample of 1382 OCs. For positions, the same conventions and reference values as Gaia Collaboration et al. (2018d) were used: in cartesian Galactocentric coordinates, the Sun is located at X = −8.34 kpc, Y = 0 pc and Z = 27 pc from the centre of the Galaxy. Similarly, we set the azimuthal velocity V c at the solar radius at 240 km s −1 . The velocity of the Sun with respect to the local standard of rest (LSR) is set to (U, V, W) = (11.1, 12.24, 7.25) km s −1 (Schönrich et al. 2010).
To calculate the uncertainties in the velocity space we used a Monte Carlo sampling of the astrometric and RV measurements and their uncertainties which we assumed Gaussian. In the following, we considered the standard deviations of the Monte Carlo samples as the uncertainties. We then ended up for the full sample of OCs with median standard deviations in the cylindrical velocities of (δV r , δV φ , δV z ) = (2.8, 3.2, 1.6) km s −1 .
For distant OCs, the astrometric uncertainties translate into large Galactic velocity errors. There are 129 OCs with standard Article number, page 5 of 16 A&A proofs: manuscript no. main deviations on one of the velocity components higher than 100 km s −1 , all of them at distances from the Sun larger than 2.2 kpc. The most extreme cases are NGC 3105 and SAI 109 which have a relative parallax error of respectively 82 and 99% and are both located at 6.9 kpc. For these two clusters, the velocity uncertainties reach several thousand km s −1 in V φ and V z . Also, as mentioned in the previous section, some OCs do not have a fully reliable mean RV because of few members, lack of information about their stability or a large dispersion among the members.
We therefore defined a High Quality Sample (HQS) composed of the OCs with a reliable mean RV (uncertainty lower than 3 km s −1 based on at least three members as defined in the previous section) and having standard deviations on (V r , V φ , V z ) below 10 km s −1 . This leaves us with a HQS composed of 418 OCs with median uncertainties on (V r , V φ , V z ) of (1.2, 1.3, 1.0) km s −1 . This is very similar to the median uncertainty of the field star sample selected by Gaia Collaboration et al. (2018d) who quoted (δV r , δV φ , δV z ) = (1.2, 1.3, 1.0) km s −1 for stars with /σ > 5. Out of these 418 clusters, an homogeneous estimation of their ages is provided by  for 411 of them. Figure 4 shows the histogram of the heliocentric distance and logarithmic ages of our different samples. Our subsample of 1382 OCs with RV measurements is 90% complete with respect to the total sample of 2017 OCs from Cantat-Gaudin et al.
(2020) up to 860 pc while our HQS is 90 % complete with respect with the total sample up to 500pc.
On the bottom panel, we can see that except for the very young clusters, all ages are represented in the HQS which indicate that our selection does not introduce a noticeable bias in terms of the age distribution for clusters older than 10 Myr. The kinematical selection based on the RV uncertainty removed most of the very young clusters, the RV of young stars being less reliable.

Kinematics of OCs compared with field stars
In this section, we compare the kinematics of the field stars sample from Gaia Collaboration et al. (2018d) with that of the OC sample. In Fig. 5 we compare the (V φ , V r ) distribution of the nearby OCs within 500 pc and 200 pc, with the field stars from Gaia DR2 closer than 200 pc and having a relative error in parallax < 0.05. In this space, field stars from the Solar neighbourhood show clear substructures in form of arches and clumps, which are associated with resonances due to non-axisymmetric features of the Galaxy, as discussed in Gaia Collaboration et al. (2018d) and subsequent works. Large clumps were found to be associated with previously known moving groups (Eggen 1958;Dehnen & Binney 1998;Chereul et al. 1999;Nordström et al. 2004), and several overdensities of field stars in this space have been linked with OCs before (Gaia Collaboration et al. 2018d). The general trend of our sample of clusters closer than 500 pc is to overlap with the higher density regions of field stars, as already seen in Soubiran et al. (2018a). The several OCs that we have gained here fall in the central part of the velocity distribution. Although we do not see a clear clumping of the clusters' distribution, it seems that OCs are most frequent in a band just between the two central arches of the field stars. Therefore, we see that OCs do not follow the exact distribution of overdensities drawn by the field stars. This was already seen with a smaller number of clusters in Gaia Collaboration et al. (2018d), where they mention that only the Pleiades and Hyades clusters are associated with overdensities, and other clusters do no show a particular overdensity in this space.
In Fig. 6 we plot the distribution of azimuthal velocities as a function of Galactocentric radius of the full sample of OCs and the HQS. Field stars in this figure were retrieved from the Gaia archive selecting stars with the same criteria as Gaia Collaboration et al. (2018d), providing a sample of 6,376,803 stars. In this space, Antoja et al. (2018) showed that the diagonal ridges can be signatures of phase-mixing after a perturbing event, or alternatively due to resonances of a barred potential or of the spiral arms. Kawata et al. (2018) argued that many of the diagonal ridges are likely related with the perturbations from the bar outer Lindblad resonance and spiral arms, which may be explained with the transient spiral arm scenario. Dias et al. (2019) and Barros et al. (2020) claimed spiral resonances are able to trap stars or OCs orbits inside the corotation radius and associated the corotation to some ridges. In our plot, the clusters seem to qualitatively follow the distribution of diagonal overdensities described by the ridges, at least the most prominent ones. The recent N-body simulations by Khanna et al. (2019) suggest that the ridges are more prominent for stars close to the plane with solar metallicity, conditions generally accomplished by a sample of OCs. Probably due to the shape of these ridges, we see a decreasing velocity gradient of V φ towards large radii, more remarkable in the plot where all the sample of clusters is shown. This is consistent with what is seen in the top middle panel of Fig. 9, where outer clusters (yellow) have smaller velocities (and in general have older ages, see the histogram). This fact, as discussed more in detail in the next subsection, produces a dip in the rotation curve traced by our sample of OCs.

Rotation curve of the Milky Way
We represent in Fig. 7 the azimuthal velocities of the full sample and the HQS OCs as a function of the Galactocentric radius, superimposing a theoretical rotation curve extracted from the axisymmetric gravitational potential described in Sect. 4. Our sample of OC data describes with good precision the rotation curve of the Milky Way for a mostly young population in the range R GC ∼[6.5,10.5] kpc. In the bottom panels, we represent the median value of the rotation velocity in bins of 400 pc, approximately. Uncertainties are computed as √ π /2 · σ / √ N. We only plot the points computed using more than 20 OCs, which we consider reliable estimations and not dominated by outliers.
Both the HQS and the full sample show similar behaviour. We see a general overlap with what is expected from the axisymmetric curve ("MWPotential2014"). Several bins in both samples depart from the theoretical curve showing a small dip towards the inner Galaxy (R GC ∼7 kpc), and a significant dip at R GC ∼9.7 kpc, which departs more than 2σ from the theoretical curve. In the full sample, the outer dip extends up to 11 kpc. In the study by Gaia Collaboration et al. (2018d), similar dips seem to be present in the sample of small Z (their fig. 13), but the inner one is only seen at positive azimuths. In the bottom panels, we overplot the unified rotation curve computed by Sofue (2020), which was obtained from a combination of data from different sources in the literature. Interestingly, this data shows a significant drop around 10 kpc, although with a different depth compared to our values.
There are several previous mentions in the literature to a dip in the region R GC ∼ 9.5 kpc (Reid et al. 2019;Barros et al. 2016;Sofue et al. 2009) whose nature is not fully understood. It has been tentatively explained by a ring density structure observed in the neutral H gas a bit further than the solar radius (R GC ∼10 kpc Nakanishi & Sofue 2016), which successfully reproduces the dip. Dias et al. (2019) calculated a corotation radius of 8.51 ± 0.64 kpc and found a similar dip in their rotation curve which can therefore be attributed to the effect of corotation, clusters located outside the corotation radius having a slower azimuthal velocity than clusters located inside. The depth of the dip is different depending on the data used, the one we observe here is compatible with the recent study by Reid et al. (2019) from star-forming regions: a small decrease from the model of ∼ 5 km s −1 . Moreover, the general decreasing slope in the outer Galaxy shown by our data is also seen with the star-forming regions. Other studies from Cepheids (e.g. Mróz et al. 2019) also show a slight decrease in V φ before and after the solar radius. McGaugh (2019) built a model that reconciles the observed stellar rotation curve with the one seen by interstellar gas, taking into account the overdensities of the spiral arms. In their Fig. 4, the modelled points (green) show a clear dip inside the solar radius, compatible with the one we see in the cluster data in Fig. 7.
Comparing with Fig. 6, these dips are a direct consequence of the shape of the most prominent diagonal ridges, clearly shown by field stars in the Gaia RVS, but also by our clusters, where a remarkable decreasing velocity gradient is seen towards the outer Galaxy. This effect is analysed in Martinez-Medina et al. (2019) using a non-axisymmetric Galactic model, that shows that the duo bar-spiral arms produce diagonal-like ridges of stars of constant angular momentum. They show how these structures tend to clump the stars in the V φ − R GC plane, pulling them up at the beginning of the ridge, and down at the end. This translates into wiggles in the rotation curve of the Galaxy. We overplot in Fig. 7 the rotation curve derived from the model by Martinez-Medina et al. (2019) 1 . Interestingly, the model resembles the two dips we observe with the OC data. We attribute to the same phenomenon the dips we observe in the cluster data, thus explained by the effect of non-axisymmetric structures in the angular momentum of the stars. Figure 8 shows the distribution of the HQS in the (V r , V φ ), (V r , V Z ) and (V φ , V Z ) planes coloured by age. Several distant old OCs with extreme velocities (V r > 40 km s −1 ) have been recently discovered. Those included in the HQS are LP 930, UBC 326 and UBC 324 all located further than 1 kpc in heliocentric distance. The most extreme OCs in total velocity ( V 2 r + V 2 φ + V 2 z ) are Ruprecht 171, Haffner 5 and Berkeley 32 already known as high velocity objects for the disc.

Age dependence of Galactic velocities
The three components of the velocities of the HQS clusters as a function of age are shown in the top panel of Fig. 9. Clearly, they show a significant growing spread as age increases, which 1 we have scaled the rotation curve of the model to fit the reference frame used in our work, described in Sect. 3, using the expression: is related with the kinematic heating and will be discussed in more details in Sect 3.4. The radial and vertical velocities are centred on zero while the azimuthal velocity is centred on the azimuthal velocity of the LSR. In addition, we can see in the middle panel of the top row that the azimuthal velocity of the oldest clusters within a given age range decreases when the radius increases. This trend is only noticeable for old clusters. Indeed, due to our kinematical selection based on the RV quality and errors of the Galactic velocities, most of the young distant clusters of our sample were removed. This bias was previously highlighted by Gaia Collaboration et al. (2018d): RVs for young stars are often less reliable due to their high temperature, their rotation or due to the fact that they are closer to the Galactic plane and not observable in spectroscopy due to heavy extinction. We can also see that inner clusters (in blue) are preferentially old, with a wide range of azimuthal velocities, a dozen of them having V φ < 220 km s −1 . For the vertical velocity shown on the right panel of the top row, the increase of the dispersion highlights the fact that the clusters reaching higher V Z are in general the oldest ones.
In the bottom panels of Fig. 9 is shown the histograms of each component of the velocities in several bins of age. The youngest clusters form a quite symmetric distribution, centered near zero for the radial and vertical components (median value of V r = −0.3 km s −1 and V Z = −1.3 km s −1 ) and near the azimuthal velocity of the LSR for the azimuthal component (median value of V φ = 239.5 km s −1 ). In the radial component, as age increases, the distribution becomes less gaussian and the peak of the histogram tends to move towards larger V r , although this is dependent on the binning of the histogram and the taken age limits. Also depending on the binning, there seems to be a sign of bimodality in the 4 older age bins, although low number statistics does not allow to reach a significant conclusion. However, independently of the binning, the oldest bin has a significantly higher maximum value (around ∼ 14 km s −1 ) with respect to the two youngest ones. Also, as age increases, the distribution becomes more and more asymmetric, with a large tail towards negative V r , and a smaller tail on the positive side becoming more pop-ulated at older ages. This effect is also independent of the histogram binning and age limits. This implies that clusters moving inwards do it with a larger range of radial motions than those moving outwards. This trend is noticeable in the whole range of ages represented in our OC population, which is representative of the youngest population in the Galactic disk. The majority of old clusters of our sample have a positive V r highlighting that OCs which are moving outwards have more chance of survival than clusters moving inwards. For the azimuthal component, the histograms in bins of ages clearly show the asymmetric drift: in the two youngest bins the distribution is centred nearly on 240 km s −1 while the mode of the histogram shifts toward lower values in the older bins. In the three oldest bin, the most probable value is around 230 km s −1 . In the third and fourth age bins, we can also note a bimodality of the distribution also seen for the radial component of the velocity, the bimodality being however independent of the binning here. As age increases, the growing spread of the vertical velocity is clearly visible on the histograms, particularly in the two last age bins, where the distribution becomes asymmetric but remains centred near zero. Figure 10 shows the distribution of V φ of the HQS sample in the X-Y plane in four different age bins. The spiral arms modelled by Reid et al. (2014) are represented by the shaded structures while the updated Cygnus arm from Reid et al. (2019) is represented by the dashed shaded structure. We note that in the two youngest bins, OCs are mainly located close to the Sun near the local arm and this spatial distribution is very different from that of the original sample presented in Cantat-Gaudin et al.
(2020) (their Fig. 8). This is a consequence of the bias of our sample stated above. The remaining nearby young clusters show a remarkable homogeneity in azimuthal velocities, with a typical velocity dispersion of 5.5 km s −1 around the Sun velocity. In the two oldest age bins, OCs are much more spatially and kinematically dispersed and do not follow the spiral arms any more.

Age velocity relation for open clusters
The increasing trend of the velocity dispersions for field stars with age is known as the age-velocity relation (AVR) and has been discussed for decades (e.g. Wielen 1977;Freeman 1987;Nordström et al. 2004). The AVR is reasonably described as a power law with some saturation at the oldest ages (Binney & Tremaine 2008;Aumer & Binney 2009;Soubiran et al. 2008), altough there is a debate about the shape of the relation (Martig et al. 2014). The AVR involves different mechanisms, such as the radial migration and heating induced by giant molecular clouds, the spiral arms (Mackereth et al. 2019) or by the bar (Grand et al. 2016), which affect differently the velocity components. The AVR study is complicated by the mixture of different stellar populations present in the solar neighbourhood, such as the thin and thick discs, which have different scale heights and different star formation and heating histories (Yu & Liu 2018;Mackereth et al. 2019). Similarly, inner and outer populations are differently affected by the dynamical effect of the bar, the spiral arms and accretion events. In addition, uncertainties in the determination of stellar ages strongly affect the AVR (Aumer et al. 2016). Ages of OCs are more reliable than ages of individual stars and determined homogeneously in . This allowed us to investigate the AVR of OCs in a large age range up to ∼2.5 Gyr (too few clusters older than this limit) and with a significant number of objects below 1 Gyr. The spatial extension is mostly confined within 200 pc from the Galactic plane, so the OC population is representative of the thin disc. It is thus interesting to compare the velocity ellipsoids of OCs and field stars in the age range where they overlap. Figure 11 shows the velocity dispersions measured in seven age bins for the 418 OCs from the HQS with an age determination. In each age bins, we fitted a Gaussian to each component of the velocities. Each point represented in Fig. 11 stands for these gaussian's standard deviations and their uncertainties. As expected, the dispersion of each component increases with age. The values and ratios are given in Table 3. The OC kinematics is characterised by a clear anisotropy in the 3 components, at all ages. The radial dispersion σ R is always significantly larger than σ φ which is significantly larger than σ z , as also seen in local field stars (e.g. Binney & Tremaine 2008;Kuijken & Tremaine 1994;Anguiano et al. 2018). The dispersion ratios remain globally stable in the different age bins, with some fluctuations within the error bars. This implies that the velocity ellipsoid keeps the same shape at all ages. Following previous studies of field stars, we fitted the increase of the velocity dispersion with age (τ) in the form of a power-law σ V ∝ τ β as suggested in Sect. 8.4 of Binney & Tremaine (2008) and by Jenkins (1992). We used a Maximum Likelihood (ML) estimator on the individual clusters of the HQS (ages younger than 2.5 Gyr), assuming Gaussian velocity errors. The results of the fits are shown in Fig 11. The grey lines represent the uncertainties on the fits: we show 100 fits taken from the final sample of the ML. The black solid line shows the best fit power law obtained with a ML, this best fit is defined as the median value of the 16000 fits performed. We found β R = 0.25 +0.05 −0.03 , β φ = 0.23 +0.03 −0.03 , β z = 0.19 +0.03 −0.03 . These values show that the heating rate of OCs is about the same in all directions, compatible with the velocity ellipsoid keeping the same shape at all ages as pointed earlier. In the following, we compare this kinematical behaviour with that of field stars. Yu & Liu (2018) measured the AVR of ∼3 500 local stars, splitting them into different subsamples depending on their Z position, and metallicity. Their metal-rich, low Z sample is comparable to our OC sample, representative of the local thin disc. For their two youngest age bins, corresponding to 1.4 ± 0.4 and 1.9 ± 0.3 Gyr respectively, they found dispersions in the different components and dispersion ratios in good agreement with what we found for OCs of similar age (see their Table 1). They fitted the heating parameters in the age range 1 to 8 Gyr and determined β R = 0.28 ± 0.08, β φ = 0.30 ± 0.09, β z = 0.54 ± 0.13, in agreement with our values in R and φ, but not in Z. Mackereth et al. (2019) measured the heating parameters in R and Z with field stars of age between 1 and 9 Gyr. They found β z = 0.50 for stars comparable to OCs in metallicity, higher than our value. They found β R to vary from 0.15 to 0.4 depending on the mean orbital radii. Recently Sharma et al. (2020) assembled several large stellar samples tracing different populations from complementary large surveys and using different age estimators. They found consistent relations for all with β R = 0.251 ± 0.006, β z = 0.441 ± 0.007, so again in agreement with our value for R but not Z. Mackereth et al. (2019) and Sharma et al. (2020) highlighted a complex relation between σ R /σ z with age, depending on metallicity, angular momentum, and height above the plane.
These comparisons of the heating parameter β (see Table 4 for more clarity) tend to show that the OC population has a dynamical evolution similar to the field stars in the radial and azimuthal directions, but not in the vertical direction. OCs seem to have a smaller heating rate in Z than field stars. The main dif-  Table 3. Dispersions of velocity components and ratios in the same bins as in Fig. 11 age interval (Gyr) N σ R σ φ σ z σ R /σ φ σ R /σ z σ φ /σ z age< 0.  Sharma et al. (2020) is the age range which extends to young ages for OCs (half of the OCs are younger than 360 million years) with few objects older than 3 Gyr, while the stellar samples typically range between 1 -10 Gyr. Although their age distributions overlap around 1 -2 Gyr, the fit of the β parameter is not performed on the same age range for OCs and field stars. It is therefore interesting that the heating rate is found similar in the Galactic plane but not perpendicular to the plane. This could mean that clusters do not reach high altitude and older ages because they are disrupted before, introducing therefore a bias in our sample. The small β z could also reflect that giant molecular clouds which are the main cause of the vertical scattering of field stars (Lacey 1984;Jenkins & Binney 1990a) are not as efficient to scatter OCs, or that the effect of the giant molecular clouds is to disrupt the OCs.
The heating of the Galactic disc and the destruction of clusters have been simulated among others by Gustafsson et al. (2016). They found that the fraction of massive old OCs, scattered into orbits with |Z| > 400 pc, is typically 0.5%. In the full  However, after our quality cuts, only 6 of such clusters remain in the HQS, making the statistics too poor to reach a conclusion on their origin.

Actions and orbital parameters
In this section we used the full 6D coordinates of the samples of clusters, to compute orbits and action variables, and analyze them as a function of the age. We used the python package galpy (Bovy 2015) for Galactic Dynamics to integrate the orbits and compute the action-angle variables. We used the axisymmetric potential MWPotential2014 implemented in galpy, which was derived by Bovy (2015) Table 3, the oldest age bin having an upper limit of 2.5 Gyr. The x value of the blue dots is the median value of the corresponding age bin. The black line shows the best fit power law obtained by using a ML on each cluster of the HQS younger than 2.5 Gyr, and the grey lines represent the uncertainties of the fits.
composed of a bulge, a Miyamoto-Nagai disc and a dark matter halo modelled by a NFW Potential. For the sake of comparison, we also used the axisymmetric potential model from McMillan (2017) to confirm our results. This model has also been fitted to the mass distribution of the Milky Way and is constituted of several components representing the cold gas disc close to the Galactic disc, both the thin and thick discs, a bulge and a dark matter halo. The detailed parameters of each component of these potentials are listed in Bovy (2015) and McMillan (2017).

Orbits
We integrated each cluster orbit with an integration step of 0.01 Myr years up to 500 Myr. We did not attempt to integrate more time (up to each cluster age) because the reliability of the results decreases significantly with time due to inaccuracies in the timedependence of the potential, amplification of the uncertainties in distance and motions, among other effects (Gaia Collaboration et al. 2018c). We computed uncertainties integrating each orbit 1000 times with a Monte Carlo sampling in the same way as described in Sect. 3. Figure 12 shows the orbits of two OCs, as an example of the results we can get using galpy. We extract the orbital parameters of the 1315 OCs for which we have an age estimation to investigate their relation with age. We represent in the left panel of Fig. 13 the evolution of the maximum altitude above the Galactic plane (Z max ) as a function of the age of the full sample of clusters and the HQS. In all the panels, the maximum height of clusters younger than 300 Myr remains constrained close to the Galactic plane. This is shown more clearly on the bottom panel where the running median of both samples increases only for an age higher than 1 Gyr. For both samples of OCs, the median age of the subsample of clusters which are reaching an altitude higher than 400 pc is greater than 1.5 Gyr. Also noticeable in the left panels of Fig. 13 is the increase of the dispersion of the maximum height of the OCs above the plane for ages older than 1 Gyr. This is usually attributed to the vertical heating of the disc: clusters are preferentially formed in the thin disc and then giant molecular clouds and spiral arms tend to scatter them away from the mid-plane (Spitzer & Schwarzschild 1951;Jenkins & Binney 1990b). This effect is consistent with what is seen in more detail in the right panels of Fig. 9, already commented in Sect. 3.
The eccentricity as a function of age in the right panel of Fig.  13 shows that : -Clusters younger than 30 Myr (log(Age) ∼ 7.5) show very low eccentricity. This is best seen in the HQS sample where the maximum value in this range is 0.018, while in the full sample there are some outliers with values up to 0.15. Very young clusters have therefore nearly circular orbits. -For clusters older than 30 Myr the dispersion of the eccentricities at a given age is large for both samples. The running median of the HQS shows an increase which is quite smooth at least up to 300 Myr (log(Age) ∼ 8.5). The full sample exhibits a similar behaviour to that of the HQS but with slightly larger mean values of eccentricities. -For ages older than 1 Gyr the dispersion in eccentricities starts to be very large for the full sample. For the HQS there seems to be a stabilization of the eccentricities around a mean value of 0.08.
This shows that OCs are born on nearly circular orbits and as their age increases, they are more prone to suffer gravitational perturbations from non-axisymmetric components. We do not find very young clusters with large eccentricities, but we find old clusters with both large and low eccentricities. We do not see a preferential location in the Galaxy, or differential characteristics between large and low eccentricity old clusters.

Action angle variables
The action-angle variables are a set of canonical coordinates which are proved to be useful to study the substructure of stars in the 6D phase space. As extensively discussed by several authors (e.g. McMillan & Binney 2008), orbital actions (J R , J φ , J z ) are integrals of motion in an axisymmetric potential, but they also provide information of non-axisymmetric perturbations. In addition, they are independent of time and are therefore more reliable to describe the orbital parameters by removing their time dependence. They have been used to describe stellar components in our Galaxy, in particular the moving groups seen in the solar neighbourhood in the (U, V) plane (Trick et al. 2019).
In an axisymmetric potential, action variables can easily be interpreted as physical quantities. The radial action J R can be used as a proxy for orbit's eccentricity or as a measure of the oscillations around the guiding radius of the object. The azimuthal action J φ is equal to the angular momentum in the vertical direction L Z , which indicates the quantity of rotation of the object around the centre of the Galaxy. Similarly, the last coordinate, the vertical action J z , can be used as a measure of the oscillations of the object around the plane of the Galaxy and therefore is a proxy for the maximum height of the object along its orbit.
All of these quantities are conserved in an axisymmetric potential but are affected by non-axisymmetric structures such as a bar, spiral arms or by a merger. We refer to Sect. 3 of Binney & Tremaine (2008) for a mathematical and comprehensive description of these variables and their meanings.
We made the computation for the 411 clusters from the HQS with known age. For comparison purposes, we made the same computation with the sample of stars within d < 200 pc in Gaia DR2 RVS, i.e. the same selection as the one done by Trick et al. (2019), which counts ∼ 350, 000 stars.
In Fig. 14 we show the distribution of radial action J R with respect to the vertical component of the angular momentum L Z , of field stars (in grey) compared with the sample of HQS clusters (coloured by the normalized density of points). In all subplots, the sample of field stars is the same as we don't know their ages, and we indicated with coloured ellipses the approximate location of the known moving groups present in the Solar neighbourhood (Antoja et al. 2008;Gaia Collaboration et al. 2018d), as analyzed in the action space by Trick et al. (2019). In this space, the location of a star is interpreted in relation to its orbital characteristics as following: (i) the "V" shape is due to the cut in distance (200 pc) made for the sample selection, (ii) stars with circular orbits are placed at (L Z , √ J R ) ∼ (1, 0), while more eccentric orbits appear at larger J R , (iii) stars close to their apocenter (pericenter) are placed in the left (right) edge of the "V".
We restrict the volume of analyzed clusters using their Galactocentric radius instead of their heliocentric distances, because of the reduced amount of clusters at d < 200 pc. Even though this made the comparison between cluster and field more difficult, the precision in the distances and velocities of the clusters is much better than for individual field stars, so we expected that the kinematic information will not be blurred by this, as it happens usually for individual stars. We made the cuts in R GC = 8.3 ± [0.2, 0.3, 0.5] pc, showed in each row. Using the information of cluster ages we were able to add an additional dimension to the figure, so we dissected the sample in four age bins [< 30,30 − 150,150 − 500,> 500] Myr, showed as columns. The aforementioned differences in the volume selection of the clusters and the field stars are visible in the top row where some of the clusters stand out of the left/right edges, which would correspond to clusters towards the Galactic centre/anticentre. From Fig. 14 we see that the bulk of clusters tends to be concentrated towards the position (1, 0), not surprising since most OCs have cold kinematics, and tend to have nearly circular orbits. It is visible in all three distance limits (rows) that there are few clusters at large eccentricities, but increasing in number at older ages. As a consequence, the density peak (marked in yellow) moves upwards towards older ages. This is particularly clear in the top row because of the large number of points.
Regarding the relation of the clusters and the moving groups, we see that there are no clusters populating the two Hercules substreams, which are stars with too eccentric orbits and low angular momentum. On the contrary, the location of the Hyades/Pleiades/Coma Berenices/Sirius moving groups resembles more the distribution of the clusters, even at larger distances (top row). We highlight that for the age bin 30 − 150 Myr (for all distance cuts) there is a significant clump of clusters towards the mid-left side (smaller L z ), populating the region where the Pleiades moving group is found. On the other hand, the Sirius moving group is populated mainly by clusters in the two older age bins. These dependencies in the population in each moving group depending on age are in the same direction as the findings by Antoja et al. (2008): in their Fig. 14, Hyades-Pleiades-Coma gather most of their young stars ( 300 Myr) and at 300 Myr there is a bump of stars related to Sirius. Hercules is associated with stars of ∼ 2 Gyr, and the distributions of all moving groups also exhibit a bump at this age, however, we almost do not find clusters at this old age, so there is no indication of this in our sample.

Summary
Thanks to the combination of Gaia and ground-based surveys and catalogues, and with new memberships from Cantat-Gaudin et al. (2020), we assembled the largest catalogue of RV velocities for OCs in order to study their kinematics. As a by-product of our study, that includes the Gaia DR2 RVS, Gaia-ESO survey, APOGEE, RAVE, GALAH and smaller catalogues, we compared the RVs from the different sources to each other to assess their typical precision and zero-point. We found RV zero points to be consistent at a level better than 1 km s −1 . The scatter of the comparisons indicate that the real precision of each catalogue is compatible with the individual uncertainties listed in it. All non-Gaia RV measurements were corrected to align them on the Gaia RVS zero point. The weighted mean RV of each star and each cluster resulted in 1 382 OCs having a RV, 38% with a highly reliable RV based on more than 3 stars and with an uncertainty lower than 3 km s −1 .
We computed both heliocentric and Galactocentric cartesian and cylindrical velocities for this sample of OCs and defined a High Quality Sample composed of 418 OCs with the most reliable velocities out of which 411 OCs have an age determination. We found that most OCs fall in a band in between the two main arches drawn by field star in the V r − V φ plane, while they seem to follow the overdensities described by the diagonal ridges in the R GC − V φ plane. The rotation curve drawn by our OCs shows two significant dips: at R GC ∼ 7 kpc, and a more prominent one around R GC ∼ 9.7 kpc. The locations and depths of these dips are in agreement with the perturbations we would expect from the non-axisymmetric components of the disc which also draw the ridges observed in the R GC − V φ plane.
With the ages of almost all the clusters from our sample, we investigated in details the age velocity relation for OCs which shows a clear anisotropy between the 3 components of the velocities. Compared with field stars studies, the heating parameter β of OCs was found to be similar in the radial and azimuthal directions, but significantly lower in the vertical direction. This low heating rate in the Z coordinate can be due to the disruption of old clusters which are the most likely to reach high altitudes above the disc or to a less efficient heating of OCs by giant molecular clouds. Although we are aware that the quality cuts we applied discarded distant clusters resulting in a bias of our sample.
We use the 6D + age information of our sample of OCs to compute and investigate orbits and action variables. We analyze the dependencies of the recovered orbital parameters as a function of age. We see that most of the clusters reach a maximum altitude above the plane during their orbits smaller than 400 pc, and only those older than 1 Gyr are able to move considerably away from the midplane, but typically less than 1 kpc. Clusters younger than 30 Myr show a very low eccentricity (∼ 0.018), and for clusters older than that, especially those older than 100 Myr, the eccentricity shows an increasing relation with age. These results show that OCs are born in circular orbits, and as age increases, they are more prone to suffer perturbations of their orbits. This is also seen after the computation of action variables, where, as age increases, the distribution in the (L Z , √ J R ) plane tends to spread beyond ∼ (1, 0). We relate our cluster distribution in this action space with the location of the known moving groups, as a function of age. We conclude that the Pleiades-Hyades-Coma moving groups seem to be more populated by young clusters, while the Sirius region seems to have a clump of clusters of age 300 Myr. No clusters are populating the two Hercules streams.
Y. Tarricq et al.: 3D kinematics and age distribution of the Open Cluster population Fig. 14. Radial action vs angular momentum (J R , L z ) distribution of field stars closer than 200 pc in Gaia DR2/RVS sample (grey), and clusters (coloured dots by the density of points). Each column shows clusters in different age bins, and each row includes clusters selected inside different ranges of Galactocentric radius, instead of heliocentric distances. We indicate the approximate location of the known moving groups analyzed by Trick et al. (2019): Hercules (green and blue), Pleiades (dark violet), Hyades (magenta), Coma Berenices (red), Sirius (yellow). We highlight the values of the three clusters Pleiades, Hyades and Coma Berenices with a star of colour violet (second column of the bottom row), pink and red (both in the last column of the bottom row), respectively, in their age-R GC panels.