Gaia EDR3 view on Galactic globular clusters

We use the data from Gaia Early Data Release 3 (EDR3) to study the kinematic properties of Milky Way globular clusters. We measure the mean parallaxes and proper motions (PM) for 170 clusters, determine the PM dispersion profiles for more than 100 clusters, uncover rotation signatures in more than 20 objects, and find evidence for radial or tangential PM anisotropy in a dozen richest clusters. At the same time, we use the selection of cluster members to explore the reliability and limitations of the Gaia catalogue itself. We find that the formal uncertainties on parallax and PM are underestimated by 10-20% in dense central regions even for stars that pass numerous quality filters. We explore the the spatial covariance function of systematic errors, and determine a lower limit on the uncertainty of average parallaxes and PM at the level 0.01 mas and 0.025 mas/yr, respectively. Finally, a comparison of mean parallaxes of clusters with distances from various literature sources suggests that the parallaxes (after applying the zero-point correction suggested by Lindegren et al. 2021) are overestimated by 0.01+-0.003 mas. Despite these caveats, the quality of Gaia astrometry has been significantly improved in EDR3 and provides valuable insights into the properties of star clusters.


INTRODUCTION
The most recent data release (EDR3) from the Gaia mission (Gaia Collaboration 2021) does not provide new data products, but instead improves upon the previous DR2 in various aspects related to the photometric and astrometric catalogues. In particular, the statistical uncertainties on parallaxes and proper motions (PM) µ have been reduced by a factor of two on average, and the systematic uncertainties are reduced even further. Already after DR2, it became possible to measure the mean parallaxes (Chen et al. 2018;Shao & Li 2019) and PM of almost all Milky Way globular clusters (Gaia Collaboration 2018; Baumgardt et al. 2019;Vasiliev 2019b) and even to study the internal kinematics of many of these systems: sky-plane rotation (Bianchini et al. 2018;Vasiliev 2019c;Sollima et al. 2019) and PM dispersion and anisotropy (Jindal et al. 2019). The improved data quality in EDR3 prompted us to reanalyze these properties, and at the same time to explore the fidelity and limitations of EDR3 itself. Since all stars in a given globular cluster have the same true parallax (with negligible spread) and share the same kinematic properties, we may use these datasets (amounting to tens of thousands stars in the richest clusters) to test the reliability of measurements in the Gaia catalogue and their formal uncertainties.
We apply a number of strict quality filters to select stars E-mail: eugvas@lpi.ru that are believed to have reliable astrometric measurements, and use these "clean" subsets to determine the properties of the cluster and the foreground populations, and individual membership probabilities for each star, in a mixture modelling procedure detailed in Section 2. We then use the selection of high-probability members to assess the statistical and systematic uncertainties on and µ in Section 3. After calibrating the recipes for adjusting these uncertainties, we proceed to the analysis of mean PM and parallaxes of clusters in Section 4. Namely, we compare the parallaxes with literature distances and find an empirical parallax offset of 0.01 mas (on top of the zero-point correction applied to each star according to Lindegren et al. 2021b) and an additional scatter at the same level. In Section 5 we explore the orbital properties of the entire cluster population, including some of the poorly studied objects with no prior measurements. Then in Section 6 we analyze the internal kinematics of star clusters: PM dispersions, anisotropy, and rotation signatures. Section 7 wraps up.

METHOD
We follow the mixture modelling approach to simultaneously determine the cluster membership probability for each star and to infer its properties, in particular, the mean parallax, proper motion, its dispersion, and other structural parameters. Mixture models have an advantage over more traditional cleaning procedures such as iterative n-sigma clipping, allowing a statistically rigorous estimation of parameters of the distribution even in cases of low contrast between the cluster and the field populations. Our procedure consists of several steps: • First we retrieve all sources from the Gaia archive that are located within a given radius from the cluster centre (this radius is adjusted individually for each cluster and is typically at least a few times larger than its half-light radius; in some cases we increase it even further to obtain a sufficient number of field stars, which are a necessary ingredient in the mixture modelling). At this stage, we do not impose any cuts on sources, other than the requirement to have 5or 6-parameter astrometric solutions (5p and 6p for short).
• Then we determine a 'clean' subset of sources that have more reliable astrometry, following the recommendations of Fabricius et al. (2021), but with tighter limits on some parameters. This clean subset excludes sources that have (a) G < 13, or (b) RUWE > 1.15, or (c) astrometric excess noise sig > 2, or (d) ipd gof harmonic amplitude > exp 0.18 (G − 33) , or (e) ipd frac multi peak > 2, or (f) visibility periods used < 10, or (g) phot bp rp excess factor exceeding the colourdependent mean trend described by equation 2 and Table 2 of Riello et al. (2021) by more than 3 times the magnitudedependent scatter given by equation 18 of that paper, or (h) flagged as a duplicated source. If the number of 5p sources satisfying these criteria exceeds 200, we use only these, otherwise take both 5p and 6p sources. The severity of these quality cuts depends on the density of stars: in the central 1 − 2 arcmin of some clusters there are virtually no stars satisfying all these criteria. For a few clusters, we slightly relax these thresholds, since otherwise they would have retained too few stars. In addition, we scale the observational uncertainties on and µ by a density-dependent factor discussed in the next section.
• Then we run a simplified first pass of mixture modelling in the 3d astrometric space only (parallax and two PM components µα ≡ (dα/dt) cos δ, µ δ ≡ dδ/dt). The distribution of all sources is represented by two (if the total number of stars is below 200) or three Gaussian components, and we use the Extreme Deconvolution (XD) approach (Bovy et al. 2011) to determine the parameters (mean values and covariance matrices) of these components, which, after being convolved with observational uncertainties (using the full 3 × 3 covariance matrix provided in the catalogue), best describe the actual distribution of sources. One of these component (usually the narrowest, except NGC 104 and NGC 362, which sit on top of the Small Magellanic Cloud with a lower PM dispersion than the cluster stars) is identified with the cluster, and the remaining one or two components -with field (usually foreground) stars.
• Then we run a full mixture model, which differs from the XD model in several aspects. We use the angular distance from the cluster centre R as an additional property of each star (besides and µ), and fit the cluster surface density profile Σ(R) by a simple Plummer model with the scale radius being a free parameter updated in the course of the MCMC run. Moreover, the intrinsic parallax dispersion of cluster members is set to zero, and the intrinsic PM dispersion is represented by a cubic spline in radius with 2−5 nodes (depending on the number of stars), with amplitudes that are also varied during the fit. Finally, we allow for spatial variation of the mean PM of cluster stars relative to the mean PM of the entire cluster. We transform the celestial coordinates α, δ, corresponding PM components, and their uncertainty covariance matrix into Cartesian coordinates and velocities on the tangent plane, using the standard orthographic projection (e.g., equation 2 in Gaia Collaboration 2018), and then convert the PM into the rotational (tangential) and radial components. The tangential motion is parametrized by a fixed profile, µt(R) = µrot 2(R/R0)/ 1 + (R/R0) 2 with R0 equal to the scale length of the Plummer density profile and a free amplitude µrot. The radial PM component of cluster members is assumed to be caused by perspective expansion (an assumption that is tested a posteriori in Section 6) and hence does not add free parameters. We use the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al. 2013) to explore the parameter space, starting from the astrometric parameters determined by XD and reasonable initial values for the remaining parameters.
• After the MCMC runs converged, we determine the membership probability for each star (including those that did not pass the initial quality cutoffs), averaging the results of classification scheme from 100 realizations of model parameters drawn from the MCMC chain. The colour-magnitude diagram (CMD) of cluster members is visually inspected to verify the outcome of the mixture model, which did not use the photometric information. The final number of members for each cluster ranges from a few up to more than 10 5 stars, with typically between 10 and 50% of stars satisfying all quality filters and hence contributing to the determination of mean properties of the cluster.
• The mean parallax and PM of the cluster and their uncertainties are taken from the MCMC chain, but these values do not take into account spatially correlated systematic errors, which we include by running an additional postprocessing step, as described in Section 6. At this step, we also re-evaluate the internal PM field of cluster members (rotation profile and radial expansion/contraction rate) using a more flexible (spline) parametrization and taking into account systematic errors.
This procedure differs from the one used in Vasiliev (2019b,c) in several aspects: (a) we typically use two rather than one Gaussian component for the foreground population in the mixture model; (b) we fit the PM dispersion profile simultaneously with other parameters, rather than determining it a posteriori from the list of members; (c) we explore the parameter space with MCMC, rather than picking up the single maximum-likelihood solution, which implies propagation of uncertainties in all nuisance parameters (including membership probabilities) into the quantities of interest such as mean PM and its dispersion. The analysis of the entire list of clusters takes a few dozen CPU hours. Similar mixture modelling approaches have been used by Sollima (2020), who focused on the outer parts of clusters and additionally used photometric information in membership classification, and Vitral (2021), who employed a fat-tailed Pearson distribution for the field star PM, since it better describes the actual distribution than a single Gaussian (for this reason, we used two Gaussian components for the field population when possible).

ERROR ANALYSIS
Gaia astrometry has been significantly improved in EDR3, both in terms of lower statistical uncertainties and better calibration. In addition, a larger number of quality criteria are available to filter out stars with possibly unreliable astrometry. Nevertheless, there are still some remaining issues with the statistical and systematic errors, which we explore in this section. These are grouped into several categories: (1) consistency between averaged values among stars of different magnitudes and spatial regions; (2) underestimated statistical errors; and (3) spatially correlated systematic errors. In these experiments, we use the samples of stars in several of the largest clusters classified as high-confidence members by the mixture model (with ≥ 90% probability, although for most stars it actually exceeds 99%).
The very first question that is natural to ask is whether the average values of parallax and PM are consistent between different subsets of stars split by magnitude or spatial location. Figure 1 shows the results of this analysis for the two largest clusters with ∼ 50 000 members: NGC 104 (47 Tuc) and NGC 5139 (ω Cen). We split the stars into three magnitude ranges, three radial intervals and four quadrants in each interval. The statistical uncertainties on the mean parallax of stars in each bin are small enough to highlight that the difference between these values exceeds these uncertainties, and is at the level ∼ 0.005 − 0.01 mas, depending on the cluster. Figure 2 shows the variation of mean parallax with magnitude for two dozen clusters (without further splitting stars into different spatial regions). Different clusters display a variety of behaviours, ranging from no systematic variation (NGC 2808, NGC 6838) to upward (NGC 5272,NGC 6341) or downward (NGC 6362, NGC 6656) trends with magnitude, or even more complex non-monotonic patterns. This diversity indicates that the parallax zero-point correction suggested by Lindegren et al. (2021b) adequately compensates these variations on average, but not necessarily for each specific region on the sky. The variations of mean parallax in different magnitude bins and the residual scatter of the overall mean parallax (deconvolved from statistical uncertainties of individual stars) are both at the level 0 − 0.02 mas depending on the cluster. These variations may not be compensated by the zero-point correction approach, since it relies primarily on sources that are either too sparse and faint (quasars) or concentrated in one region of the sky (stars in the Large Magellanic Cloud -LMC), so should be considered as an unavoidable additional systematic uncertainty. We now turn to the analysis of its spatial variation.
It is known that the Gaia astrometry contains spatially correlated systematic errors associated with the scanning law, which vary mostly on the scale 0.5 • , but have some correlations over larger angular scales as well. Although the prominence of these 'checkerboard patterns' is significantly reduced in EDR3 compared to DR2 (see, e.g., Figure 14 in Lindegren et al. 2021a), they are not completely eliminated. The lower panels of Figure 1 show the 2d maps of mean parallaxes in NGC 104 and NGC 5139, which indeed fluctuate on a scale of 10 − 20 . In the first approximation, these spatially correlated errors in astrometric quantities χ ≡ { , µ} can be described by covariance functions Vχ(θ) ≡ χi χj, which depend on the angular separation θ between two points i, j (measured in degrees). Lindegren   . Spatial covariance function for parallax as a function of separation between sources V (θ), estimated in four richest clusters in different ranges of magnitudes: from G = 13 up to G = 18 (left), G = 19 (centre) or G = 21 (right panel), note that each next range includes the previous one. Pairs of stars are binned into 7 − 15 distance bins; the uncertainties in each bin are driven by measurement errors, which increase with magnitude. In all panels, the values of V (θ = 0) are significantly different from zero, but they increase with the inclusion of fainter stars.
functions in a number of ways, primarily based on the sample of ∼ 10 6 quasars, which are found across the entire sky (except the regions close to the Galactic disc), but are relatively faint. They find V (θ) 140 µas 2 on scales ∼ 1 • , and possibly a few times higher in the limit of small separation, although with a large statistical uncertainty limited by the number of close pairs of quasars. On the other hand, LMC stars are brighter and more dense, allowing them to extend the covariance function to smaller separations, which turns out to be ∼ 3× smaller than for quasars. Maíz Apellániz et al. (2021) carried out an independent analysis of quasars and LMC stars, but restricted their sample to quasars brighter than G = 19. They also obtained lower values for V in the range 50 − 80 µas 2 for separations smaller than a few degrees, and their LMC sample yielded a yet lower V ∼ 10 − 50 µas 2 for θ 1 • . All these findings suggest that the spatial correlations are less prominent for brighter stars, but the dependence of the covariance function on magnitude has not yet been quantified (although Fardal et al. 2021 found that the amplitude of the 'checkerboard pattern' in DR2 increases by a factor of a few from bright to faint stars, qualitatively similar to what we see in EDR3).
We explored the spatial covariance of parallaxes and PM in the four richest clusters (NGC 104, 5139, 6397 and 6752), using only the 5p stars that pass all quality criteria, but considering different magnitude ranges. Figure 3 shows that the binned covariance functions have a similar behaviour for all clusters, and its limiting value at θ = 0 does depend on the magnitude cutoff, increasing towards fainter magnitudes. The value for bright stars (13 < G < 18) is V (0) 50 µas 2 , consistent with the findings of Maíz Apellániz et al. (2021), and a few times large if we consider all stars. Due to a finite spatial extent of cluster members, the estimated covariance function drops to zero or negative values at separations θ 0.2 − 0.3 • comparable to the size of the cluster, which does not reflect its true behaviour on large scales (in other words, the systematic error may have spatial correlations on larger scales, but they would be the same for all cluster members and hence simply shift the mean , which is subtracted before computing V ).
We thus augment our estimates of V (θ) with the ones from Lindegren et al. (2021a) and Maíz Apellániz et al. (2021) for θ 0.2 • , while introducing an overall magnitude-dependent normalization factor. Our approximation is where G is the magnitude of the brighter star in the pair. We do not use the same functional form as Maíz Apellániz et al. (2021) since the latter is unphysical (a valid covariance function must have a non-negative angular power spectrum, i.e., Legendre integral transform), but instead reproduce the overall trends with a different function.
To conduct a similar analysis for the PM covariance function, we relied on the assumption (tested later in Section 6) that the radial component of PM is entirely caused by perspective effects, thus after subtraction of these it should be zero on average. Unfortunately (in this context), unlike parallax, the intrinsic dispersion of PM is significantly larger than measurement errors, thus reducing the statistical significance of the estimated Vµ(θ). Nevertheless, we obtain a reasonably consistent limiting value Vµ(0) 400 [µas yr −1 ] 2 for all clusters and magnitude ranges; combined with the larger-scale trends from Lindegren et al. (2021a), we approximate These covariance functions yields a systematic uncertainty 0.011 mas and µ 0.026 mas yr −1 in the limit of small separations, which should be viewed as the irreducible systematic floor on the precision of parallax and PM measurement for any compact stellar system. The overall uncertainty (statistical and systematic combined) for a given selection of cluster members is derived using the method described in Vasiliev (2019c), and may be smaller than these values if a cluster spans more than a fraction of a degree on the sky. Maíz Apellániz et al. (2021) discuss a similar method for combining statistical and systematic uncertainties, but in their equations 5-7, the overall error is dominated by the largest values of V (θ), whereas our method essentially sums up the values of the inverse covariance matrix, and therefore puts more emphasis on stars with smallest spatial covariances. In practice, though, the difference should be minor.
Apart from confirming the systematic error, we also find that the formal statistical uncertainties are somewhat underestimated. This is most easily manifested in the distribution of parallaxes for cluster members: the intrinsic spread of is negligibly small even for the closest clusters ( 10 −3 mas), and hence we should expect the uncertainty-normalized deviations from the mean parallax for each star, ( i − )/ ,i, to follow a standard normal distribution. Figure 4 shows the distribution of these measurements for one of the richest clusters, NGC 5139 (ω Cen), split into several groups by stellar magnitude and distance from the cluster centre. It is clear that even for the "clean" subset of stars that satisfy all quality criteria, the distribution has more prominent tails than a Gaussian, at least in the central regions with high density of stars. We may quantify this effect by assuming that the actual (externally calibrated) measurement uncertainty˜ is given by the scaled formal uncertainty, summed in quadrature with  Histograms are split by distance range (from top to bottom row) and magnitude (from left to right, with the last column showing 6p sources, which are typically fainter than G = 20). If the formal uncertainties were correctly estimated, the histograms should have followed a standard normal distribution (shown by gray dotted lines), however in practice, the distribution is broader by a factor η 1.05 − 1.2, depending on the density of sources (and hence on radius).  Figure 5. Parallax error inflation factor η as a function of source density and magnitude, for a dozen clusters with ≥ 5000 members. For each cluster, we use only the clean subset of stars (5p in the first three panels, and 6p in the last one) split into several equal-number bins in radius; the horizontal axis shows the mean density of stars in each bin, and the vertical axis -the width of the best-fit Gaussian fitted to the distribution of uncertainty-normalized deviations from the mean parallax for stars in the given bin ( i − )/ ,i , similar to the one shown in Figure 4 for one particular cluster. The scaling factor η that the uncertainties should be multiplied by to match the standard normal distribution is generally larger in higher-density central regions; the gray dotted lines show the trend η = (1 + Σ/Σ 0 ) ζ with Σ 0 = 10 stars arcmin −2 and ζ = 0.04. The first and the last panel exceed this trend line on average, and a better fit could be obtained by adding in quadrature a systematic error sys = 0.01 mas (first panel -bright stars) or lowering the value of Σ 0 (last panel -6p sources). . PM dispersion profiles for NGC 3201 (left and centre panels) and NGC 5272 (M 3, right panel) as a function of distance, split into several magnitude bins. Left panel: when using nominal measurement uncertainties (the average µ for each bin is printed in the legend), the deconvolved intrinsic dispersion is higher for fainter stars, indicating that their uncertainties are underestimated. In the remaining two panels, we use the same prescription for the density-dependent error inflation factor η as for the parallax, which produces profiles that are in reasonable agreement between different magnitude bins. Shaded areas show the 68% confidence intervals on σµ(R) from the MCMC runs of the full mixture model pipeline, in which this scaling factor was also applied and faint stars have been excluded. 2 ,ext = η 2 2 + 2 ,sys .
(3) Figure 5 shows the error inflation factor η estimated in several magnitude bins from a dozen clusters, as a function of source density. Despite some scatter, the overall trends are consistent between all clusters, and indicate that the formal uncertainties should be scaled by η 1.1 − 1.15 in regions with high source density, but η 1 in lowerdensity regions. For brighter stars, η is somewhat higher, but one can compensate for this by an additive systematic error ,sys 0 − 0.02 mas depending on the cluster (the additive errors for individual clusters are shown on Figure 2). These values are comparable with those reported by Fabricius et al. (2021, their figures 19, 21), although they did not study the variation with source density and did not apply all quality criteria that we used here, thus their η are somewhat higher on average. El-Badry et al. (2021) and Zinn (2021) find a similarly mild ( 1.2) error inflation factor from analysis of parallaxes of wide binaries and astroseismically calibrated stars, respectively (after applying some quality cuts similar to ours). We also examined the performance of the composite quality filter ("astrometric fidelity") suggested by Rybizki et al. (2021), but found that it fails to remove many astrometrically unreliable sources in dense central regions, producing much broader distributions of normalized parallax deviations with η ∼ 1.5 − 2.
The right panel of Figure 5 shows that the scaling factor η is slightly higher for 6p sources that pass all other quality filters. We may approximately describe its dependence on source density Σ as η = (1 + Σ/Σ0) ζ ; the values of Σ0, ζ and sys for various subsets are reported in Table 1.
It is natural to expect that the PM uncertainties could be underestimated by a similar factor, but this is more difficult to test empirically, since the intrinsic PM dispersion is non-negligible for most clusters with a sufficiently high number of stars, and is indeed among the free parameters in the mixture model. One possibility is to check whether the errordeconvolved PM dispersion is the same when computed from different magnitude ranges. Figure 6 demonstrates that when using formal uncertainties from the Gaia catalogue without any correcting factors, the internal PM dispersion appears to be higher for fainter stars, indicating that their uncertainties are likely underestimated. On the other hand, when adopting the same prescription for the PM uncertainty scaling factor η as for the parallax, the inferred dispersion is usually consistent between different magnitude ranges. Based on these experiments, we adopted this scaling prescription for the entire mixture model fitting procedure. We also checked that the results are largely insensitive ( 0.01 mas yr −1 difference) to 5% variations in the adopted value of η 1.15 in the highest-density regions (and corresponding changes at lower densities).
Nevertheless, to avoid possible biases caused by clusterto-cluster variations of the scaling factor, we conservatively used only those stars from the clean subset that have sufficiently small uncertainties (usually these are the brighter ones). Namely, we select only stars with µ < κσµ/(η − η0), where κ = 0.2 is the tolerance parameter, η(Σ) ≥ 1 is the default density-dependent error scaling factor, and we set η0 = 0.9. The idea is that in high-density central regions, η is not only higher, but also more uncertain, and thus we reduce the maximum acceptable statistical error of stars in this region to limit the bias from incorrect assumptions about this scaling factor. Obviously, this filter already requires the knowledge of the intrinsic dispersion σµ, so we apply it iteratively, using the profiles σµ(R) obtained in previous runs. All stars in the clean subset that do not pass this filter are still used in the astrometric fit, but are convolved with a fixed (previous) σµ(R) profile instead of the one inferred during the fit, so do not bias its properties. If the number of stars with small enough uncertainties is too low ( 50), we do not attempt to infer the PM dispersion and instead assume a fixed profile guided by the line-of-sight velocity dispersion, or simply zero for distant or low-mass clusters. A similar cut on statistical uncertainties on the HST -derived PM was applied by Watkins et al. (2015), who removed all stars with uncertainties larger that 0.5σµ (this is a more stringent criterion than we use, but their PM uncertainties are typically at the level 0.05 mas yr −1 , which is lower than most Gaia stars).

MEAN PARALLAXES AND PROPER MOTIONS
Having established the necessary adjustments and correction factors for the mixture modelling procedure, we now discuss its outcomes pertaining to the global properties of the clusters -mean parallaxes and PM. Figure 7 shows the comparison of Gaia parallaxes with the distances compiled by Baumgardt & Vasiliev (2021) from various literature sources. In this plot, the parallaxes and their uncertainties are computed by simple error-weighted averages over member stars (with individual zero-point corrections for each star), without accounting for correlated systematic errors. As expected, the majority of points lie near the line D = 1, but the deviations from this line are often significantly larger than could be expected from the formal statistical uncertainties alone. Moreover, there is a general tendency of parallax to be slightly larger than 1/D, illustrated by a large fraction of clusters with 5 < D < 20 lying above the D = 1 line. This may indicate that the parallax zero-point correction suggested by Lindegren et al. (2021b) is slightly overshooting, although it is also possible that there are some general systematic biases in the literature distances. Before exploring this question in detail, we discuss a few specific anomalies.
A number of clusters at low Galactic latitudes have systematically lower parallaxes than implied by the literature distances, but they are in highly-reddenend regions and it is plausible that the CMD-derived distances are biased. In particular, for Pal 6 and Pal 7 (IC 1276), it appears that using the parallax distance and slightly adjusting the extinction coefficient produces a better fit to the CMD than the literature values. Parallaxes of NGC 4147 and Terzan 1 are also significantly smaller than 1/D, and despite rather large uncertainties, the difference is statistically unlikely, possibly indicative of undetected problematic sources in Gaia. Finally, for some bright, low-extinction clusters such as NGC 6809 (M 55), NGC 6266 (M 62), NGC 288 and NGC 7089 (M 2), the disagreement between CMD-derived and parallax distances exceeds the statistical uncertainty by a factor of few, and cannot be explained other than by the systematic variation of parallax zero-point across the sky.
To explore the possible systematic biases and additional uncertainties, we performed the following analysis. We assume that the deviations of the measured distance moduli 1 µi and parallaxes i from their true values for i-th cluster follow normal distributions: where µ,i are formal uncertainties on the distance moduli 1 they are traditionally denoted by the same symbol µ as the PM.
from the literature, ,i are the Gaia statistical uncertainties, µ,sys, ,sys are the additional global systematic uncertainties, and ∆µ, ∆ are global offsets (same for all clusters). For each choice of these four global parameters, we compute the joint likelihood of measured values µi, i by marginalizing over the unknown true distance modulus: We then maximize the joint likelihood of all clusters in a Monte Carlo simulation, restricting the sample to clusters with ≥ 100 stars and low extinction, E(B − V ) < 0.5 (we also tried a number of other quality cuts, but the results were qualitatively similar). Figure 8 shows the posterior distribution of these four parameters, and the gray-shaded region in Figure 7 shows the samples from this distribution on the parallax vs. distance relation. The offset and additional error in distance modulus are consistent with zero, though the uncertainties are at the level of a few×10 −2 mag. On the other hand, the offset in Gaia parallax is ∆ 0.01 ± 0.003 (that is, zero-point corrected parallaxes are slightly too high), and the additional systematic error in parallax is ∼ 0.01 -a very similar value to the one obtained in the previous section from two completely unrelated arguments (spatial covariance function and the variation of mean parallaxes in different magnitude bins). This value likely represents the true limit of the parallax precision in Gaia EDR3.
The overcorrection of the parallax zero-point by the Lindegren et al. 2021b prescription has been noted in several other studies: for bright stars (G < 13), Riess et al. (2021) and Zinn (2021) find that the Gaia parallaxes are overestimated by ∼ 0.015 mas, while Huang et al. (2021) find an overestimation of ∼ 0.009 mas for stars fainter than G = 14, with a further position-dependent variation of ∼ 0.01 mas. Our estimated offset is qualitatively similar to these studies. Only in special circumstances it may be possible to estimate the zero-point offset more precisely at the location of a given cluster: Chen et al. (2018) combined the sparse quasars with the much more numerous stars of the Small Magellanic Cloud (SMC) to derive corrected parallaxes for NGC 104 (47 Tuc) and NGC 362 from Gaia DR2. Their values 104 = 0.225 ± 0.007, 362 = 0.117 ± 0.007 agree well with our measurements (0.232 and 0.114 respectively). Recently Soltis et al. (2021) derived the parallax for NGC 5139 (ω Cen) from Gaia EDR3 to be 5139 = 0.191±0.004, which, unsurprisingly, agrees with our value 0.193, as well as with a number of other distance estimates from the literature. However, we note that despite a good agreement, this value still carries a systematic uncertainty 0.01 mas, corresponding to a 5% distance uncertainty, significantly larger than they optimistically assumed, and reducing the precision of their calibration of cosmological distance indicators. We still need to wait until further Gaia data releases to bring down the systematic uncertainty to competitive levels.
Comparing our parallax values with the ones derived by Shao & Li (2019) from Gaia DR2, we find a general agreement within error bars, after correcting for the mean parallax   ). The vertical axis shows the Gaia parallax multiplied by the distance: values above 1 correspond to the parallax distance being smaller than the literature value, and vice versa. The vertical error bars take into account the statistical uncertainties on both the mean parallax and the distance, but the horizontal error bars for the distance are not displayed; all clusters with D < 0.2 are shown. Parallaxes of individual stars that contribute to the mean parallax of a cluster are corrected for the colour-and magnitude-dependent zero-point offset following the recipe from (Lindegren et al. 2021b). Points are coloured according to the number of cluster members, and empty points indicate clusters with higher reddening E(B − V ) ≥ 0.5, for which many photometric distance determination methods may be less reliable, or with fewer than 100 stars. The agreement between Gaia parallaxes and other distance measurements is fairly good, but at large D the parallaxes seem to be higher on average than 1/D, and the scatter is larger than the statistical uncertainties. The gray-shaded region shows Monte Carlo samples from a statistical relation (4) reproducing this offset and scatter, with parameters shown in Figure 8. offset in DR2 ∆ −0.03 (i.e., DR2 parallaxes are smaller on average). Only a few clusters showed a statistically significant disagreement, e.g., NGC 5272 (M 3) ( our − S19 0.085), NGC 6544 (0.07), NGC 7099 (M 30) (0.07). The parallax values for the few brightest clusters considered in Maíz Apellániz et al. (2021) are in complete agreement with our measurements, being derived from the same EDR3 catalogue. Likewise, the mean PM of clusters in EDR3 are in good agreement with the ones derived by Gaia Collaboration (2018), Baumgardt et al. (2019) and Vasiliev (2019b) from Gaia DR2 data: for 80% clusters, the total difference in both PM components is within 0.1 mas yr −1 , comparable to the systematic uncertainty of DR2 (0.066 mas yr −1 per component). Clusters with the largest PM difference usually have very few stars or are located in dense regions: AM 4 (∼ 1.7 mas yr −1 ), FSR 1735 (∼ 1.0), UKS 1, Terzan 6 (0.8), Terzan 5, Djorg 1 (0.6). Recently Vitral (2021) independently derived mean PM for ∼ 100 globular clusters from Gaia EDR3, which are in a very good agreement with our results.
The previous analysis involved the mean parallaxes derived using the statistical uncertainties alone. Given the several independent pieces of evidence for the additional systematic error, in the subsequent analysis we use the parallaxes and PM computed with full account for spatially correlated systematic errors, following the method of Vasiliev (2019c); these values are reported in Table A1. The resulting uncertainties on parallax and PM are usually at the level dictated by the covariance functions (1, 2), that is, 0.011 mas, µ 0.026 mas yr −1 per component; however, they could be slightly smaller for clusters with a large spatial extent, or larger for clusters with too few members. The mean values computed with and without accounting for systematic errors are usually very close, however, there are a few exceptions (mainly for large clusters). For instance, in NGC 104 (47 Tuc) the relatively sparsely populated outer regions have coherently higher parallaxes (lower left panel in Figure 1). Since the covariance function drops with distance, the pairs of stars in opposite 'corners' of the cluster have a larger contribution to the overall mean value computed while accounting for spatial correlations. On the other hand, the simple average parallax weighted purely with statistical uncertainties of individual stars is dominated by the more numerous stars in the inner part of the cluster, which happen to have somewhat lower parallaxes than the outskirts. As a result, the simple statistical average parallax is lower than the value computed using spatial correlations by ∼ 0.005 mas (note the offset be-  (4) between Gaia parallaxes and literature distances, sampled from a Monte Carlo chain. µ and ∆µ are the additional distance modulus error and its offset, both are consistent with zero. is the systematic uncertainty on Gaia parallaxes that needs to be added to random uncertainties to make them statistically consistent with distance moduli from the literature, while ∆ is the parallax offset (positive values indicate that Gaia parallaxes are on average higher than 1/D).
tween the red and the gray lines and points in the upper left panel of Figure 2). However, these offsets in NGC 104 and a few similar clusters are well within the overall systematic uncertainty on the mean parallax or PM.
The newly derived mean PM are dominated by systematic uncertainties at the level ∼ 0.025 mas yr −1 for most clusters (unless the number of members is below ∼ 100 or the contrast between the cluster and the field stars in the PM space is low, in which case statistical uncertainties may be higher). However, we are ultimately interested not in the value of the PM, but of the transverse physical velocity v ⊥ ≡ µ D, and the uncertainty in distance usually is the dominant limiting factor in the precision of the velocity, as illustrated in Figure 9. Overall, the velocity uncertainty is of order ∼ 2 − 20 km s −1 for the majority of clusters, except the most distant ones.
Out of 157 objects in the Harris (1996Harris ( , 2010 catalogue, we could not determine the PM for only a few clusters which are located in highly extincted regions and are not visible to Gaia: 2MASS-GC01, 2MASS-GC02, GLIMPSE01 and GLIMPSE02. We have also added a number of recently discovered globular cluster candidates to our list: FSR 1716 (Minniti et al. 2017 . Uncertainty in the transverse velocity v ⊥ ≡ µ D as a function of heliocentric distance D has two components: distance uncertainty µ D is shown by circles (individual estimates where available, or assuming a rather optimistic 5% relative error otherwise), and PM uncertainty µ D is shown by crosses (taking into account systematic errors µ 0.025 √ 2 mas yr −1 , which usually dominate over statistical errors). For most clusters, the first factor is more important. Colours show the number of cluster members with reliable astrometry for each object.
1 (Muñoz et al. 2012), BLISS 1 (Mau et al. 2019), bringing the total count to 170. Many of these additional objects are poorly studied and lack line-of-sight velocity measurements, and the nature of them is not well established. We shall see below that a few of these are located at low Galactic latitudes and move within the disc plane, so may well be open rather than globular clusters.

ORBITS OF GLOBULAR CLUSTERS
Given the full 6d phase-space coordinates of clusters, one may examine their orbital properties or the distribution in the space of integrals of motion: this has become a popular exercise especially after the advent of precise PM from Gaia DR2 (e.g., Binney & Wong 2017;Myeong et al. 2018;Massari et al. 2019;Piatti 2019;Forbes 2020;Kruijssen et al. 2020;Pérez-Villegas et al. 2020;Bajkova et al. 2020). The entire population of clusters can be used to probe the gravitational potential of the Milky Way, using Jeans equations or distribution function-based dynamical models (e.g., Watkins et al. 2019;Posti & Helmi 2019;Vasiliev 2019b;Eadie & Juric 2019). However, the outer parts of the Galaxy, where this exercise is most useful, are subject to the non-equilibrium distortions caused by the recent passage of the Large Magellanic Cloud (LMC), which perturbs the velocities by as much as few tens km s −1 and can bias the inference on the gravitational potential (Erkal et al. 2020a,b;Cunningham et al. 2020;Petersen & Peñarrubia 2021). Although the effect of the LMC can be accounted for in dynamical models, as illustrated by Deason et al. (2021) in the context of a simple scale-free model for halo stars, we leave a proper treatment of globular cluster dynamics for a future study. Figure 10 shows the orbital properties of the entire population of globular clusters: semimajor axis, eccentricity and Horizontal axis shows the semimajor axis, vertical -eccentricity, and colour -inclination (prograde in red, retrograde in blue, and polar orbits in green). Each cluster is shown by a cloud of points sampled from the measurement uncertainties, which in most cases are dominated by distance uncertainties. We note that the changes in orbit parameters in different Galactic potentials sometimes exceed the measurement uncertainties, especially in the outer part of the Galaxy.
inclination, computed in two variants of static Milky Way potentials: McMillan (2017) or Bovy (2015) (these properties remain qualitatively similar if we use other reasonable choices for the potential). Each cluster is shown by a cloud of points representing the measurement uncertainties of its 6d phase-space coordinates, which, as discussed above, are mostly dominated by distance uncertainties (also propagated into the transverse velocity). For the 9 clusters lacking line-ofsight velocity measurements, the missing velocity was drawn from a global distribution function, resulting in a plausible mean value with a large uncertainty 100 km s −1 . An alternative depiction is the rhomboid action-space diagram (e.g., Figure 5 in Vasiliev 2019b), which shows the same information in a different projection of the 3d space of integrals of motion. Note however that peri-and apocentres (equivalently, eccentricity and semimajor axis) are computed by numerical orbit integration, whereas actions are computed in the Stäckel approximation; in practice, this distinction is unimportant.
Objects located in the same region of the plot and having similar colours are close in the 3d integral space and may be physically related, such as the population of globular clusters associated with the Sgr stream. Note that the proximity in the integral space does not imply that the objects line up on the same path on the sky plane, therefore the association with the stream should be examined in the space of observables -celestial coordinates, distances, PM and lineof-sight velocities. This connection was explored in a number of papers, e.g., Law & Majewski (2010b), Bellazzini et al. (2020), Arakelyan et al. (2020); however, these studies all relied on the old model of the Sgr stream from Law & Majewski (2010a), which was conceived before the more recent observations of its trailing arm, and does not match its features (primarily the distance to the apocentre). When using the model of the Sgr stream from Vasiliev et al. (2021), which adequately matches all currently available observational constraints, we find that only the following clusters can be unambiguously associated with the Sgr stream: four clusters in the Sgr remnant -NGC 6715 (M 54, which sits at its centre), Terzan 7, Terzan 8, Arp 2, three clusters in the trailing arm -Pal 12, Whiting 1 and NGC 2419, and one faint cluster Ko 1, which previously lacked PM measurements (and still has no line-of-sight velocity data), but now coincides with the stream in position, distance and both PM components. Interestingly, it sits in the region where the second wrap of the trailing arm intersects with the second wrap of the leading arm, so is consistent with both interpretations. The line-of-sight velocity is expected to be rather different: 100 to 150 km s −1 for the trailing arm or −100 to 0 km s −1 for the leading arm. Ko 1 was previously conjectured to belong to Sgr stream by Paust et al. (2014) together with its sibling Ko 2; however, the latter does not match the stream neither in distance nor in PM. A few distant outer halo clusters -Pal 3, Pal 4 and Crater -might also be associated with the Sgr debris scattered beyond the apocentre of the trailing arm, but they do not match the model track in at least one dimension (though the model might be unreliable beyond 100 kpc, not being constrained by observational data).
Other conspicuous features in Figure 10 include the populations of high-eccentricity clusters with semimajor axes in the range 6 − 20 kpc, which are believed to be associated with an ancient merger of a satellite galaxy on a nearly radial orbit (Myeong et al. 2018). Finally, a significant fraction of clusters in the inner part of the Galaxy (with semimajor axes below 6 − 8 kpc) have low eccentricity, prograde disc-like orbits (coloured dark red/purple on the plot). Some of the recently discovered objects in the outer part of the Galaxy also share these characteristics. We computed the reflex-corrected PM (assuming that the distance is known to sufficient accuracy), and for several clusters at low Galactic latitudes (|b| < 25 • ), the PM component perpendicular to the Galactic plane (µ corr b ) is significantly smaller than the parallel component (µ corr l ). We may conclude that the orbits of these clusters necessarily stay close to the disc plane, even if the line-of-sight velocity is not known (these cases are marked by italic). This subset of disc-like clusters includes BH 140, BH 176, BLISS 1, ESO 93-8, Ko 2, Mercer 5, Pfleiderer 2, Segue 3. Some of these objects may well be old open clusters rather than globular clusters. The orbit of Ryu 059 also lies close to the disc plane but is retrograde and highly eccentric (although its line-of-sight velocity is not known, the reflex-  Figure 11. Slope of the radial PM component ξ ≡ µ R /R (in units of mas yr −1 per degree). Left panel shows the measured values for 10 clusters with uncertainties ξ smaller than 0.05, which agree fairly well with values expected from perspective effects. Right panel shows the distribution of deviation of measured values from expectations, normalized by measurement uncertainties ξ , for ∼ 150 clusters with at least 100 stars. The red curve shows the standard normal distribution, which adequately describes the actual histogram of deviations (though the latter has a slight excess of objects around zero, for which the uncertainties ξ might be overestimated).
corrected sky-plane velocity already exceed 400 km s −1 ), with the estimated apocentre radii exceeding 100 kpc. It is quite plausible that the distance to this cluster is overestimated -a similar story happened with Djorg 1, for which a 30% downward distance revision by Ortolani et al. (2019) made its orbit much less eccentric and more realistic. Finally, Kim 3, Laevens 3, Muñoz 1, Ryu 879 have high-inclination orbits.
As is clear from Figure 10, the orbital parameters often have significant uncertainties, but these are strongly correlated, therefore quoting the confidence intervals on r peri/apo or eccentricity makes little sense without a full covariance matrix, or better, the full posterior distribution. Rather than providing these quantities in a tabular form, we provide a Python script 2 for computing orbits in any given potential, sampling from the possible range of initial conditions for each cluster; the integrations are carried out with the Agama library for galactic dynamics (Vasiliev 2019a). Figure A1 in the Appendix shows the derived PM dispersion and rotation profiles for more than 100 clusters that have a sufficiently large number of stars and are not too distant, comparing our results with other studies.

INTERNAL KINEMATICS
We explored the sky-plane rotation signatures in all clusters, using the method detailed in the Appendix A6 of Vasiliev (2019c), which takes into account spatially correlated systematic errors. Namely, we fitted a general linear model with N + 3 free parameters to the PM field of stars with high membership probability, where the free parameters are the two components of the mean PM, the slope ξ of the radial PM component (µR = ξ R), and N amplitudes of a B-spline representation of the tangential PM component µt, with N = 2 − 4 depending on the number of stars. As discussed in that paper, it is not possible to combine the analysis of spatially correlated systematic errors and probabilistic membership, so we used it as a post-processing step after the main MCMC run, and accounted for uncertain membership by considering 16 realizations of the subset of cluster members, selecting stars in proportion to their membership probability. The PM dispersion profile was kept fixed, as the previous analysis indicated that it is little affected by the spatial correlations (unlike the mean PM field).
The radial PM component is expected to be caused entirely due to perspective contraction or expansion: for a cluster at a distance D moving with line-of-sight velocity v los , the expected µR at an angular distance R (measured in degrees) is ξ R, with ξ expected = −v los /D × (π/180 • /4.74). Figure 11, left panel, compares the measured and expected values of ξ for 10 clusters with small measurement uncertainties ξ , demonstrating a very good agreement. The cluster with the largest amplitude of the perspective expansion is NGC 3201, for which the value of ξ is measured with 6% relative error; for NGC 5139 (ω Cen) the relative error is ∼ 12%. The right panel of the same figure shows the histogram of differences between measured and expected values, normalized by the measurement uncertainties, which closely follows the standard normal distribution (though with an excess of points around zero, for which the uncertainties ξ might be overestimated). We stress that these uncertainties take into account spatially correlated systematic errors: if we use only statistical errors, ξ would be considerably smaller for rich clusters, and deviations between measured and expected ξ would exceed (3 − 4) ξ for quite a few objects. This exercise validates our approach for treating correlated systematic errors and gives credence to the similar analysis of rotation signatures in the PM.
The uncertainty in the PM rotation profile is typically at the level 0.01−0.015 mas yr −1 for sufficiently rich clusters, being dominated by systematic errors. We detect unambiguous rotation (at more than 3σ level) in 17 clusters, with further Table 2. Rotation signatures in star clusters detected in this work (last column) and in some previous studies based on Gaia data: Gaia Collaboration (2018), Bianchini et al. (2018), Vasiliev (2019c), Sollima et al. (2019). "+" indicates a firm detection, "?" -a tentative detection, and "−" stands for no clear signature. : that study considered both PM and line-of-sight velocities, and all these clusters have inclination angles exceeding 60 • , i.e., the rotation signal is mostly seen in the line-of-sight velocity field (though we do see weak signatures in the PM field of the last two objects). We also excluded Terzan 5 due to a small number of stars satisfying our quality cuts.
Turning to the analysis of PM dispersion profiles, we compare them with the profiles derived from Gaia DR2 by Vasiliev (2019c) and Baumgardt et al. (2019), finding generally a good agreement, with the present study having somewhat smaller uncertainties due to improvements in Gaia astrometry. In addition, the inferred PM dispersion profiles can be compared with the HST -derived PM dispersions in the central parts of 22 clusters studied in Watkins et al. (2015) and 9 clusters studied in Cohen et al. (2021). Unfortunately, due to the strict quality cutoffs adopted in the present study, there is very little (if any) spatial overlap between Gaia and HST measurements from Watkins et al. (2015), but in general, the PM dispersion profiles agree remarkably well. Some discrepancy is seen in NGC 5139 (ω Cen), where our PM dispersion profile is lower (although there are almost no stars in our clean subset to anchor it in the crowded central parts), NGC 6341 (M 92), where Gaia σµ is slightly higher, NGC 6535, where both profiles have large uncertainties but Gaia is lower, and NGC 6681 (M 70), where the Gaia σµ profile is significantly lower. The last case is the most puzzling discrepancy, which appears to be robust against variations in the assumed error inflation factor η or various quality filters. The 9 clusters with HST measurements from Cohen et al. (2021) have a larger spatial extent and agree very well with Gaia in all cases except NGC 6355 and NGC 6401, where Gaia σµ is slightly higher, and NGC 6380, where Gaia dispersion is ∼ 15% lower.
The PM dispersion profiles can also be compared with lineof-sight velocity dispersion profiles from various spectroscopic studies; here we use three such datasets - Kamann et al. (2018) used MUSE IFU in central regions of some clusters, while Ferraro et al. (2018) and Baumgardt et al. (2019) provide wide-field coverage. The conversion from σµ to σLOS involves distance, and is one of the methods for its determination. Using the average distances from the literature, we generally get good agreement between PM and line-ofsight dispersions, but in some cases the discrepancy is significant and may indicate some problems in either dataset, or alternatively, calls for a revision of the distance. For wellpopulated clusters, the discrepancies are usually less than 10%. The most notable outliers are NGC 1904, NGC 5272 (M 3), NGC 6388, where Gaia PM is ∼ 10 − 15% higher for the adopted distances; NGC 6304, NGC 6553, NGC 6626, NGC 6779 (M 56), where Gaia PM is ∼ 15% lower; and NGC 6681 (M 70), which is problematic as described above.
Although we used isotropic PM dispersion in the mixture model by default, in the richest clusters we were able to explore the anisotropy by fitting the radial and tangential PM dispersions separately. Figure 12 shows the radial profiles of PM dispersion anisotropy, σt/σR − 1, for 16 clusters with sufficient number of stars and satisfying the consistency checks that the PM dispersion and anisotropy agree between stars in different magnitude ranges. There is a considerable diversity in the anisotropy profiles: half of these clusters have predominantly radially anisotropic PM dispersion, a few others have tangential anisotropy, and NGC 5139 (ω Cen) transitions from being radially anisotropic in the inner part to tangentially anisotropic in the outskirts. The difference between σR and σt is typically at the level 10−20%, except NGC 7078 (M 15) -a core-collapsed cluster with a rather strong radial anisotropy in the outer part. Reassuringly, the line-of-sight velocity dispersion outside the core of M 15 matches well the tangential component of PM dispersion (since both mainly reflect the tangential component of 3d velocity), while the radial component of PM dispersion is noticeably higher (another reason why the line-of-sight velocity should never be called "radial velocity"!).
PM anisotropy was explored by Watkins et al. (2015) in the central parts of 22 clusters observed by HST, who found deviations from isotropy at a few per cent level. However, as discussed above, there is almost no spatial overlap between HST and the clean Gaia subset. Our findings can be more directly compared to Jindal et al. (2019), who explored PM anisotropy in 10 clusters using Gaia DR2. Our results agree for all clusters except NGC 6656 (M 22), which they found to be radially anisotropic but we do not see a strong evidence for this, and NGC 6397, which appeared to be isotropic in their analysis but weakly radial in our study. The radial anisotropy in NGC 3201 has also been detected by Bianchini et al. (2019) using Gaia DR2, consistent with our measurement.

SUMMARY AND DISCUSSION
Gaia EDR3 is a quantitative rather than qualitative improvement upon the revolutionary DR2 catalogue, yet its precision is materially better. In the first part of this study, we scrutinized the fidelity and accuracy of Gaia astrometry with a variety of methods, using the data from a dozen richest globular clusters, after applying a number of stringent filters to remove possibly unreliable sources.
• Formal statistical uncertainties on parallax and PM adequately describe the actual errors in regions with low stellar density, but appear to be underestimated by 10 − 20% in higher-density regions; Table 1 provides the suggested correction factors.
• Spatially correlated systematic errors in parallax and PM are considerably lower than in DR2: for bright stars, the residual systematic error in is at the level 0.01 mas (a fourfold improvement), while for stars fainter than G = 18 it may be twice higher. The systematic uncertainty in PM is ∼ 0.025 mas yr −1 , or 2.5× better than in DR2.
In the second part, we re-examined the kinematics of almost all Milky Way globular clusters, derived the mean parallaxes and PM, galactic orbits, and analyzed the internal rotation, dispersion and anisotropy profiles of sufficiently rich clusters. While the improvements in precision for mean PM with respect to studies based on DR2 (Gaia Collaboration 2018, Baumgardt et al. 2019, Vasiliev 2019b) is substantial, the precision of the phase-space coordinates is typically limited by distance rather than PM uncertainties. Gaia parallaxes will eventually deliver the most precise distance estimates for nearby clusters, but at present, the systematic uncertainty at the level 0.01 mas limits their usefulness in practice. On the other hand, the internal kinematics (PM dispersions and rotation signatures) have also improved considerably, and agree well with various independent estimates. We find evidence of rotation in more than 20 clusters, and measure the PM dispersion profiles in more than a hundred systems, down to the level 0.05 mas yr −1 , i.e. at least 2× better than in DR2. These data can be used to improve dynamical models of clusters and provide independent distance estimates, which are examined in Baumgardt & Vasiliev (2021).

ACKNOWLEDGEMENTS
EV acknowledges support from STFC via the Consolidated grant to the Institute of Astronomy. We thank A.Riess and G.Gontcharov for valuable discussions. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/ dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

DATA AVAILABILITY
We provide the catalogues of all Gaia sources in the field of each cluster with their membership probabilities and the tables of radial profiles of PM dispersion and rotation amplitudes, available at https://zenodo.org/record/ 4549397, as well as the summary table of mean parallaxes and PM and scripts for computing cluster orbits in a given potential, available at https://github.com/ GalacticDynamics-Oxford/GaiaTools.