The Hubble Space Telescope UV Legacy Survey of Galactic Globular Clusters. XXIII. Proper-motion catalogs and internal kinematics

A number of studies based on data collected by the $\textit{Hubble Space Telescope}$ ($\textit{HST}$) GO-13297 program"HST Legacy Survey of Galactic Globular Clusters: Shedding UV Light on Their Populations and Formation"have investigated the photometric properties of a large sample of Galactic globular clusters and revolutionized our understanding of their stellar populations. In this paper, we expand previous studies by focusing our attention on the stellar clusters' internal kinematics. We computed proper motions for stars in 56 globular and one open clusters by combining the GO-13297 images with archival $\textit{HST}$ data. The astro-photometric catalogs released with this paper represent the most complete and homogeneous collection of proper motions of stars in the cores of stellar clusters to date, and expand the information provided by the current (and future) $\textit{Gaia}$ data releases to much fainter stars and into the crowded central regions. We also census the general kinematic properties of stellar clusters by computing the velocity-dispersion and anisotropy radial profiles of their bright members. We study the dependence on concentration and relaxation time, and derive dynamical distances. Finally, we present an in-depth kinematic analysis of the globular cluster NGC 5904.


INTRODUCTION
Over the past 20 yr, photometric and spectroscopic data have radically changed our picture of Galactic globular clusters (GCs). One of the most baffling discoveries is the presence of multiple stellar populations (mPOPs) in GCs (see Anderson 1997;Lee et al. 1999;Bedin et al. 2004;Carretta et al. 2009a,b;Piotto et al. 2015;Renzini et al. 2015;Bastian & Lardo 2018;Gratton et al. 2012Gratton et al. , 2019Cassisi & Salaris 2020, and references therein). To shed light on the formation and evolution of mPOPs, the "HST Legacy Survey of Galactic Globular Clusters: Shedding UV Light on Their Populations and Formation" (GO-13297, PI: Piotto;Piotto et al. 2015) was devised. UltraViolet (UV) and optical data of this and other Hubble Space Telescope (HST ) programs have allowed the creation of color-magnitude (CMDs) and color-(pseudo-)color diagrams, which, in turn, have provided essential elements for better understanding the mPOP phenomenon in Galactic GCs (e.g., Monelli et al. 2013;Milone et al. 2017).
Other recent important observational findings concern the internal kinematics of GCs. Studies on the internal motions of these systems have revealed that GCs are characterized by complex internal kinematic properties including velocity anisotropy, rotation, and partial energy equipartition (see the review of Varri et al. 2018, and references therein). These results have provided the motivation for new theoretical studies aimed at building the theoretical framework necessary to interpret these observational findings (e.g., Bianchini et al. 2017Bianchini et al. , 2018bTiongco et al. 2016aTiongco et al. , 2017Tiongco et al. , 2018Tiongco et al. , 2019Breen et al. 2017Breen et al. , 2021Szölgyén et al. 2019;Pavlík & Vesperini 2021, 2022. Most of what we have learned on this topic comes from internal motions in the plane of the sky obtained with HST and Gaia data. These two space observatories have different characteristics that make them best suited for specific, yet complementary, investigations of GCs. The Gaia mission has revitalized astrometry. The availability of high-precision proper motions (PMs) over the entire sky has enabled a large variety of investigations, for example, large-scale structures in the Galaxy, Galaxy kinematics, stellar streams, tidal tails (e.g., Gaia Collaboration et al. 2018a;Ibata et al. 2019Ibata et al. , 2021. The internal kinematics of GCs have benefited from the Gaia PMs as well (e.g., Bianchini et al. 2018a;Jindal et al. 2019;Evans et al. 2021), * We dedicate this paper to the memory of Prof. Ivan R. King, member of this collaboration who passed away on August 31, 2021, whose pioneering work led to a deeper understanding of globular clusters.
but these Gaia-based analyses are focused on (and limited to) bright stars outside the centermost regions. Crowding (Pancino et al. 2017) and faintness (i.e., access to low-mass stars) are two hurdles that will be complex (or impossible) to overcome even in the next Gaia data releases. However, there is important information in the cores of GCs and in their faint members that is necessary to properly characterize GCs, and one of the few ways of obtaining these necessary data is with HST.
Here we combine the wealth of information available in the HST archive to compute high-precision PMs with the goal of analyzing the internal motions within GCs (following Bellini et al. 2014). The internal kinematics have a lot to tell, not only about GCs as a whole, but also about their mPOPs. Indeed, the present-day trends in the velocity-dispersion and anisotropy radial profiles are the results of the different initial conditions of first-and second-generation stars. Thus, measuring the internal motions of the mPOPs can help us to shed light on key aspects of the mPOP formation and evolution (e.g., Mastrobuono-Battisti & Perets 2013, 2016Hénault-Brunet et al. 2015;Tiongco et al. 2016bTiongco et al. , 2019Vesperini et al. 2013Vesperini et al. , 2021Sollima 2021). This is one of the remaining goals of the GO-13297 program. We present the PM catalogs for the 56 globular and one open clusters targeted by the GO-13297 program, and an overview of the kinematics of their brightest (and more massive) members. To showcase the quality of our catalogs, we have also analyzed in great detail the internal kinematics of NGC 5904.
Together with the PM catalogs, we release several photometric catalogs that are useful for selecting highquality objects for studying the internal motions. Recently, Nardiello et al. (2018) published the final version of the photometric catalogs of this project. The photometry presented in this work is not meant to replace that published by Nardiello et al. (2018). Indeed, although the photometric precision is comparable, the completeness of the catalogs presented here is lower because they contain only sources with positions measured in at least two epochs to enable the determination of the PM. Also, we only release the photometric catalogs of the images in the filters/cameras actually used for the PM computation.

DATA SETS AND REDUCTION
One of the goals of this project is to provide highprecision, homogeneous PMs from heterogeneous HST data (i.e. observations with different cameras and programs) for stars in the central regions of GCs. We aimed at computing reliable PMs for as many stars as possible, with particular care for sources in the most crowded regions of the field of view (FoV) and for stars as faint as possible that we can detect with the available HST archival images. No two data sets are alike, thus a careful data reduction has been specifically tailored to each GC. Although the data-reduction process is detailed in various papers (Bellini et al. 2018;Libralato et al. 2018aLibralato et al. , 2019Nardiello et al. 2018), we provide here a brief overview and the main differences in our reduction with respect to these works. Two clusters were processed independently as part of other publications: NGC 362 (Libralato et al. 2018a) and NGC 6352 (Libralato et al. 2019). For these, we refer the reader to the corresponding papers for detailed descriptions of the data reduction.
We made use of flt images (which are dark and bias corrected, and have been flat-fielded, but not resampled) taken with the Wide-Field Channel (WFC) and High Resolution Camera (HRC) of the Advanced Camera for Surveys (ACS), and with the Ultraviolet-VISible (UVIS) and InfraRed (IR) channels of the Wide-Field Camera 3 (WFC3) taken before 2019 1 . In the case of ACS/WFC and WFC3/UVIS data, images were pipeline-corrected for the charge-transfer-efficiency (CTE) defects, as described in Anderson & Bedin (2010). As discussed in Bellini et al. (2014), not all filters/cameras are suitable for PMs. However, some parts of the reduction process (e.g., the second-pass photometry described below) take advantage of a large number of images to better detect faint sources. For this reason, we chose to include all images at our disposal in the first part of the data reduction 2 .
Our data reduction is a combination of first-and second-pass photometric stages. First-pass photometry was used to create an initial set of positions and fluxes for the brightest and most isolated sources in each exposure via effective-point-spread-function (ePSF) fitting through a single wave of finding. The ePSFs were specifically tailored to each image, starting from the publicly-available, spatially-variable library of HST ePSFs 3 . Source positions were also corrected for geometric distortion by means of the distortion solutions provided by Anderson & King (2004, Anderson (2016), Bellini & Bedin (2009) and Bellini et al. (2011).
1 For NGC 6121, we excluded data from the GO-12911 program, which are part of a more detailed analysis in progress on this cluster. 2 The data sets used in this work are collected at 10.17909/gajx-kf45. All our data products are available at MAST as a High Level Science Product via 10.17909/jpfd-2m08. See also: https://archive.stsci.edu/hlsp/hacks. 3 https://www.stsci.edu/∼jayander/HST1PASS/LIB/ Bright, unsaturated stars in the single-image catalogs were cross-matched with those in the Gaia Data Release 2 (DR2) catalog (Gaia Collaboration et al. 2016. This step was necessary to set up a common, pixel-based reference frame with specific axis orientation (X and Y axes toward west and north, respectively) and pixel scale (40 mas pixel −1 ). The centers of the clusters (from Goldsbury et al. 2010 with the exception of NGC 5897 and NGC 6791 for which we used the coordinates from the Harris 1996 catalog, 2010 edition, and Simbad database 4 , respectively) were placed at a specific position, for example (5000, 5000), so to always have positive master-frame coordinates (the exact coordinates of the center of the clusters is provided in the header of the published catalogs). As in Bellini et al. (2018) and Libralato et al. (2019), a master frame was created for every filter, camera and epoch. Then, we iteratively cross-identified the same stars in all images, and applied six-parameter linear transformations to transform the stellar positions in each single-image catalog onto the master frame. Once on the same reference system, the positions and instrumental magnitudes (rescaled to the magnitude of the longest exposure in each epoch/camera/filter) were averaged to create a preliminary astro-photometric catalog.
The second-pass photometry was obtained with the software KS2 (Anderson et al., in preparation; see also Bellini et al. 2017a). KS2 makes use of all images at once to increase the signal of faint objects that would otherwise be undetected in a single image. Starting from the brightest sources, KS2 progressively finds fainter stars, and measures their position and flux via ePSF fitting after all detected close-by neighbors have been subtracted from the image. This step is particularly important in crowded environments like the cores of GCs. We run KS2 separately for different epochs, by grouping data taken 1-2 yrs apart, to retain stars that have moved by more than one pixel from one epoch to another (see Bellini et al. 2018). KS2 allows us to define a specific set of images that can be "stacked" together and used to find sources in the FoV. There is not a one-size-fits-all solution for selecting these reference exposures. For each cluster, we selected a combination of cameras/detectors/filters that provided a good compromise between using a large number of exposures to more easily identify faint stars, and avoiding spurious detections (caused, for example, by the mix of filters with very different wavelength coverages).
Instrumental magnitudes in the KS2-based catalogs were converted on to the VEGA-mag system. Photometry obtained with ACS/WFC F606W, ACS/WFC F814W, WFC3/UVIS F336W or WFC3/UVIS F438W filters was registered on to the VEGA-mag system by computing the zero-point difference with the corresponding photometry in the catalogs of Nardiello et al. (2018). The zero-point takes into account for the normalization to 1-s exposure time, the aperture correction and the VEGA-mag zero-point. Photometry with all other cameras/filters was calibrated as described in Bellini et al. (2017a) by means of drz (for ACS/HRC) or drc (for ACS/WFC or WFC3/UVIS) images and the official aperture corrections and VEGA-mag zeropoints 5 .
Finally, KS2 provides the position and flux of all the detected sources in the raw reference-frame system of each image (Bellini et al. 2018). We made use of these KS2-based single catalogs to compute our PMs.

PROPER MOTIONS
PMs were computed following the procedures and caveats described in Bellini et al. (2014). Briefly, positions in the KS2-based raw catalogs were corrected for geometric distortion and then transformed (with six-parameter linear transformations) onto the same reference-frame system defined in Sect. 2. Only cameras/detectors/filters best suited for astrometry 6 (see Bellini et al. 2014) were used in the PM computation. Then, the positions as a function of time were fit with a least-squares straight line. The slope of this straight line provides a direct estimate of the PM of the source. The PM errors are the uncertainties of the PM fit obtained by using the actual residuals of the points about the fitted line (Bellini et al. 2014) For each source, the six-parameter linear transformations used in the PM computation were obtained by using a set of close-by, bright, well-measured cluster 5 See the resources provided here: https://www.stsci.edu/hst/instrumentation/acs/data-analysis, https://www.stsci.edu/hst/instrumentation/wfc3/dataanalysis/photometric-calibration. 6 Besides the filters for which we do not have an ad-hoc geometricdistortion correction, we did not use filters bluer than F336W for the WFC3/UVIS detector and F330W for the ACS/HRC camera, respectively. We also excluded all WFC3/IR data given the worse resolution of the WFC3/IR detector. 7 Objects with a peculiar motion that cannot be modeled by a simple straight-line fit, like wide binaries, would result in large PM errors. However, these objects are hard to be discerned from single stars with poorly-measured proper motions from the information in our catalogs alone. A systematic search and accurate PM estimate for objects with very-peculiar motions would require an ad-hoc analysis, which is outside the scope of this project.
stars. Thus, our PMs are computed relative to the bulk motion of each GC at that given specific location in the FoV, and the cluster PM distribution is centered on the origin of the vector-point diagram (VPD). In addition, we also provide the PM zero-point needed to transform these relative PMs on to an absolute reference system (Appendix A). Another important feature of our PM derivation is that any signature of internal cluster rotation in the plane of the sky is removed from the cluster stars, but it is present, with the opposite sign, in all other sources (Bellini et al. 2017b). Thus, we cannot directly infer cluster rotation from the kinematics of the cluster members. The same argument is also valid for the parallax effect (e.g. Libralato et al. 2018a,b).
Small spatially-variable and magnitude-dependent systematic errors are present in our uncorrected PMs. As in Bellini et al. (2018), we notice two main types of systematic errors: • a low-frequency effect correlated with the temporal baseline, as well as number and types of images, used to compute the PMs. To remove this time-dependent systematic, we divided our sample in N sub-groups based on the PM temporal baseline. N varies from cluster to cluster due to the heterogeneous data sets used. Then, we computed the median PM of each sub-group, which should be zero by construction. If it was not zero, we subtracted this median PM value to the PM of each star in the sub-group; • a high-frequency, spatially-and magnitudedependent systematic error. By construction, the average PM of cluster stars should be zero regardless their magnitude and location in the field. This is not always true locally, mainly because of a combination of CTE and geometric-distortion residuals. These residual systematic errors were corrected using the median PM of the closest, both spatially and in magnitude, N well-measured cluster members (target star excluded). The closeness criterion and the number of reference stars N were tailored to each cluster to reach the best compromise between mapping the variations as locally as possible and the need for large statistics. For very-bright(faint) objects, we set up a magnitude threshold above(below) which reference stars are used for the correction instead of a fixed ∆mag. This was done to increase the statistics at the extreme ends of the magnitude range.
The errors that we report in our catalogs for the PMs thus corrected include the propagated contributions  from the uncertainties in the corrections themselves. When no enough reference stars were available for the high-frequency correction, this correction was not applied. Figure 1 shows maps of raw and a-posteriori, locallycorrected PMs for the GC NGC 5272. The FoV has been divided into square cells of 100 WFC3/UVIS pixels per side. In each cell, we selected the 50 well-measured, cluster stars closest to the center of the cell, and computed the average PM in each direction. Panels (a) to (d) present the local PM maps obtained by means of the raw, uncorrected PMs, while panels (e) to (h) show the corrected PM maps. In each row, the two leftmost panels refer to stars brighter than the F606W instrumental magnitude equal to −10 (signal-to-noise ratio of ∼100; m F606W ∼ 21.6), while the maps for fainter stars are shown in the two rightmost panels. The comparison between the top and bottom panels clearly highlights the effectiveness of the a-posteriori corrections. However, some residual high-frequency systematic errors are still present among the faintest objects. Thus, caution is advised when using these objects.
The raw PMs of some clusters, mainly those including the ACS/WFC data of the GO-14235 program (PI: Sohn), present larger systematics related to uncorrected CTE. The CTE that affects the HST detectors has worsened over time, and the official pipeline is not always able to completely correct it. An example is shown for NGC 1261 in the top panels of Figure 2, where it is clear that (i) bright and faint stars have different systematic errors in the PMs, and (ii) the local PM map of faint stars presents a discontinuity in the FoV along the chip gap of the GO-14235 data set. For these specific cases, we applied an additional a-posteriori correction prior to that for high-frequency systematics. Briefly, we divided our sample of well-measured stars into four magnitude bins. In each sub-sample, we computed the average PM (in each direction) of cluster members in 125-pixel-wide bins along the direction perpendicular to that of the CTE systematic on the local PM map. The corrections to the PM of each star were computed by interpolating among these binned values. As for the other a-posteriori corrections, the errors of this CTE-related correction are included in the corrected-PM error budget.
The steps described above generally remove the majority of systematic errors included in the raw PMs. However, these corrections are not perfect, especially for very faint stars, and we advise users to carefully check PMs for magnitude/color/spatial systematics on a cluster-by-cluster basis. Figure 3 provides an overview of the PM catalog of NGC 6652. This GC is located in the outer Bulge, pro-jected on the Baade's Window (Rossi et al. 2015), and in foreground of the Sagittarius Dwarf spheroidal. The VPD obtained with our corrected PMs is presented in panel (a). On the basis of their locations in the VPD, we arbitrarily define three groups of objects: cluster stars (points within the red circle), Sagittarius-Dwarf members (points within the blue circle), and Bulge stars (objects within the green ellipse). In the CMD in panel (b), we highlight members of NGC 6652 in black (stars within the red circle in the VPD), Bulge objects in green and stars associated with the Sagittarius Dwarf in azure. This is a simple example of one of the possible applications of our PM catalogs. Corrected PMs in each coordinate as a function of m F606W are plotted in panels (c) and (d), and as a function of (m F606W − m F814W ) color in panels (e) and (f). Only cluster members are shown. The red points (with error bars) are the median values of the PMs in 0.5-mag bins. The azure line is set to zero. These plots do not show any significant magnitude-or color-dependent systematics in our corrected PMs, but very blue and red objects hint at the presence of some residual color-dependent systematics. However, stars with these extreme colors are very faint (see panel b), and we expect both PMs and their corrections still to be affected by small systematic residuals. We stress, though, that no quality selections other than the membership were applied here. Finally, panel (g) shows the 1D corrected-PM error as a function of m F606W . The red, horizontal line is set at the median 1D PM error of bright, well-measured stars, i.e., 25.6 µas yr −1 .
A description of the final PM and photometric catalogs (together with some caveats about their usage) is provided in Appendix B, Tables A2 and A3, respectively. The precision of our PMs varies from cluster to cluster and depends on temporal baseline, as well as on the number and depth of images. Figure 4 presents pie charts of the temporal baselines and the number of images used to compute PMs in all our catalogs, respectively. About ∼33% of the PMs are computed with only 10-20 images and a temporal baseline between 7 and 9 yr, i.e., by combining only GO-10775 and GO-13297 data. The remaining PMs result from the mix of heterogeneous data sets. For bright, well-measured stars, the median raw-PM precision is between about 7 and 60 µas yr −1 , while that of corrected PMs ranges between about 13 and 120 µas yr −1 .

INTERNAL KINEMATICS OF STELLAR CLUSTERS
As a benchmark for our PM catalogs, we here provide velocity-dispersion and anisotropy radial profiles of the clusters in our project. The kinematic profiles  of the cores complement those in the outskirts derived with the Gaia Early Data Release 3 (EDR3) catalog (e.g., . We also include the kinematic profile of the open cluster NGC 6791. We started by selecting well-measured objects in each epoch/filter/camera combination as follows (see Table A3 for the explanation of the parameters): (i) magnitude rms lower than the 90-th percentile of the distribution at any given magnitude. An object with a magnitude rms better than 0.01 mag is always included in the well-measured sample, while all sources with a magnitude rms larger than 0.15 mag are excluded.
(ii) QFIT value larger than the 90-th percentile of the distribution at any given magnitude (note that the closer to 1 the QFIT parameter is, the better the PSF fit). Again, we retained all objects with a QFIT larger than 0.99 and discarded those with a value lower than 0.75.
(iii) |RADXS| value lower than the 90-th percentile of the distribution at any given magnitude, but keeping all objects with a |RADXS| lower than 0.01 and rejecting those with a value larger than 0.15; (iv) the photometric N phot u /N phot f ratio (see Table A3) is greater than 0.75; (v) o < 1; (vi) flux at least 3.5σ above the local sky.
These (empirically-derived) thresholds were chosen as a compromise between rejecting bad measurements and keeping a large sample of stars for the subsequent analyses. To avoid crowding bias, we measured the 90-thpercentile trends for magnitude rms, QFIT and RADXS in a region outside the core of each cluster where sources are more isolated and applied these cuts to all stars across the FoV. When not enough stars were available outside the core of the cluster (for example, for the photometry obtained with ACS/HRC full-frame or WFC3/UVIS subarray images), all stars were used, regardless of their location in the FoV.
For each epoch, all criteria from (i) to (vi) have to be fulfilled in at least two filters if a star has been measured in at least two filters, otherwise in the only filter through which it was been detected. Finally, an object that has passed all previous selections is defined as "wellmeasured" if it passes all criteria in at least two epochs.
The PMs also provide useful parameters for selecting trustworthy sources for the kinematic analysis. In addition to all the previous criteria, we also removed all objects that have an astrometric N PM u /N PM f ratio (see Table A2) smaller than 0.8-0.9 (the exact value changes from cluster to cluster), χ 2 µα cos δ and χ 2 µ δ larger than 1.25-1.5, PM error larger than 0.5 mas yr −1 , or for which the a-posteriori PM correction was not computed. Finally, we also excluded stars with a PM error larger than f times the local velocity dispersion σ µ of a sub-sample of well-measured, close-by cluster stars. For each cluster, we compared the velocity-dispersion radial profiles of RGB stars obtained by varying f from 0.5 to 0.9, with steps of 0.1, finding a general good agreement (within 1σ) between the inferred kinematics. If enough stars were left after all our quality selections, we used f = 0.5. For 9 GCs (NGC 5053, NGC 5466, NGC 5897, NGC 6144, NGC 6366, NGC 6496, NGC 6535, NGC 6584 and NGC 6717), we used a value of f = 0.8, while for NGC 6981 we set f = 0.9. We included these 10 clusters in our analysis, but because their PM errors are of the order of their σ µ , we advise caution in the interpretation of their velocity-dispersion profiles.
All these criteria are designed to obtain a good compromise between statistics, quality and completeness. The only two exceptions to the strategy described above are NGC 362 and NGC 6352, for which we adopted the quality selections described in Libralato et al. (2018a) and Libralato et al. (2019), respectively.
GCs are old, collisional systems that, after many two-body relaxation times, present a (partial) degree of energy equipartition. Because of the heterogeneous mass ranges covered by our PM catalogs, we restricted the analysis of the kinematic profiles to stars brighter than the MS turnoff along the sub-giant and red-giant branches (SGB and RGB, respectively). Since in the post-MS evolutionary stages the evolutionary lifetimes are quite shorter than the core H-burning one, and in any case shorter than the typical two-body relaxation time, SGB and RGB stars can be safely considered as having the same kinematic mass (see also Sect. 6.3).
The velocity-dispersion profile of each cluster was obtained by dividing the sample of massive cluster members 8 into equally-populated radial bins, with at least 100 stars per bin. In case of low statistics, we lowered the number of stars per bin to ensure at least 3 radial bins, regardless of the number of stars within each bin. When the statistics allowed, we also imposed a radial bin comprising only the centermost 5-10 arcsec. The velocity dispersion in each radial bin was computed similarly to what is described in Raso et al. (2020), i.e., by maximizing the following likelihood: where (v rad, n , v tan, n ) are the radial and tangential components of the PM of the n-th star, ( rad, n , tan, n ) are the radial and tangential components of the PM uncertainty of the n-th star, (v rad , v tan ) are the radial and tangential mean motions of the cluster, and (σ rad , σ tan ) are the radial and tangential velocity dispersions of the cluster. We also computed the combined velocity dispersion σ µ in each radial bin using the same likelihood in Eq. 1 but with σ rad = σ tan = σ µ . We used the affineinvariant Markov Chain Monte Carlo (MCMC) method emcee  to sample the parameter space and obtain the posterior probability distribution functions (PDFs) for σ µ , σ rad and σ tan . We run the MCMC chain with 20 walkers for 5000 steps, then rejected the first 200 steps. The best-fit values correspond to the medians of the PDFs, while the corresponding errors are defined as the average between the 16-th and 84-th percentile about the median. Finally, the velocity dispersions were corrected as described in Watkins et al. (2015a, but see also Appendix A of van de Ven et al. 2006) to take into account the maximumlikelihood estimators being biased and underestimating the true dispersion of the velocity distribution. The difference between the corrected and uncorrected σ µ is, on average, ∼0.6%, it never reaches 3% and it is more important for bins with fewer than 100 stars.
Figures A2-A5 in Appendix C show the result of our analysis for the 57 stellar clusters in the GO-13297 project. These profiles are also available at our website 9 . For each cluster, the m F606W versus (m F606W −m F814W ) CMD is shown on the rightmost panel. Well-measured members of the GC brighter than the MS turnoff (highlighted by the azure, dashed horizontal line) are plotted as red points, all other objects that passed the quality criteria are shown as black points.
The velocity dispersion σ µ as a function of distance from the center of the cluster is presented in the bottomleft panel. Black, filled points are obtained from this work, while black, open points refer to the measurements in the GC database of Holger Baumgardt 10 obtained with the Gaia EDR3 PMs, which were presented by . The only exception is NGC 6791, which is not included in the work of . For this cluster, we independently derived Gaia-based σ µ using our tools 11 (Table 1). The black dashed lines are set at the core (r c ; obtained using the definition in Spitzer 1987, see Eq. 1-34, which is similar to the definition of the King scale radius r 0 ) and the projected half-light (r h ) radii provided in the Baumgardt's database. For NGC 6791, we used r c and r h from Dalessandro et al. (2015) and Kamann et al. (2019), respectively.  used all stars at their disposal regardless of their magnitude. Thus, the Gaiabased velocity dispersions for close-by clusters were derived from stars with various masses and, because of the effects of energy equipartition, could be systematically higher than our HST -based profiles. However,  argued that the Gaia uncertainties for faint stars are likely underestimated, and included a scaling factor for the PM errors to obtain consistent σ µ values between different magnitude intervals. Furthermore, the authors applied various quality cuts to their samples prior to the determination of the PM velocity dispersions (see their Sect. 3), which likely excluded from the fits faint (and so low-mass) objects with large PM uncertainties. For this reason, we choose to directly compare our HST profiles in the cores with these Gaia velocity dispersions outside the cores. Overall, we find a good agreement between HST and Gaia profiles, with a few exceptions. NGC 5466 and NGC 6981 have HST PM uncertainties of the same order of the σ µ (we used a large value of f for the analysis), thus the inferred values of σ µ should be interpreted cautiously. For NGC 6934, most of the Gaia stars used by  have PM errors larger than the intrinsic σ µ of the cluster, and the Gaia measurements might not be completely reliable. The HST profile of NGC 6304 seems higher than the one from the Gaia PMs. However, the HST -based kinematics are in agreement with the lineof-sight (LOS) velocities in the Baumgardt's database.
We fit these points with a 4-th order monotonicallydecreasing polynomial. The coefficients of the 4-th-order polynomial are obtained with a maximum-likelihood approach. For each cluster in Figures A2-A5, the blue line in the bottom-left panel is obtained with the best-fit (median) values of the polynomial fit, while the cyan lines are 100 random solutions of the polynomial fit. We 11 We considered only cluster stars (selected by means of PMs and CMD) within 750 arcsec from the center of NGC 6791. We rejected all sources that had a re-normalized unit weight error (RUWE) greater than 1.4, an astrometric excess noise larger than 0.4, a number of bad along-scan observations exceeding 1.5% of the total number of along-scan observations, or a 2D PM error worse than 0.3 mas yr −1 . We also excluded all objects fainter than the MS turnoff or brighter than G = 13. used these polynomial functions to derive the σ µ at the center of each cluster, at r c and at r h . We provide all these values in Table A4. The polynomial fits show again the agreement between HST and Gaia data for most of the clusters.
There are a few points that are outliers with respect to the polynomial-fit predictions. Most of these outliers refer to the velocity dispersion in the innermost bins. A rise of the velocity dispersion in the innermost region can be a proxy for the presence of an intermediate-mass black hole (Greene et al. 2020, and references therein), but also for issues related to crowding/blending. Indeed, if two sources are blended and confused as one, or if the light contamination from the neighbors is high (which can still occour even with our data reduction), the position measured can be shifted from the real position. The offset is different for every image/filter/camera. Hence, the net result is that its PM is likely larger than what it should be, thus increasing the σ µ of the sources in the very-crowded region (see the discussion in Bellini et al. 2014). Finally, the shape and the abrupt drop of the polynomial functions in the outermost parts of the FoV do not correspond to physical effects and are just plotted for completeness.
The anisotropy (σ tan /σ rad ) as a function of distance from the center of the cluster is presented in the top-left panels. Anisotropy is discussed in , but their values are not publicly available. Thus, only HST -based data are shown in the top-left panels of Figs. A2-A5 (black, filled dots). The red, dashed, horizontal line is set to 1 and marks the isotropic case. Most clusters are isotropic in the core, but a few objects show a radial anisotropy outside about one r h . We will discuss the kinematic anisotropy in detail in the next Section. Jindal et al. (2019) analyzed 10 GCs using Gaia DR2 data and made publicly available their velocity dispersion and anisotropy radial profiles. Although their profiles do not cover the centermost region observed by our HST data, we find an agreement at the 1σ level between their σ µ and anisotropy values for the nine clusters in common with our data set.
We also compared our kinematic profiles with those of Watkins et al. (2015a), see Fig. 5, where a previous version of the HST PM catalogs was used (Bellini et al. 2014) 12 . There is an overall agreement between the profiles at the 3σ level. The profiles of Watkins et al. (2015a) for NGC 1851, NGC 6441 and NGC 7078 are generally higher than those in our work. For NGC 7078, we also do not see the drop of the velocity dispersion in the centermost region as seen by Watkins et al. (2015a). A similar discrepancy is found when comparing the profiles of NGC 2808, NGC 6681 and NGC 6715, although at to lesser extent. The reason for these discrepancies might be the better treatment of crowding in our data reduction and/or the additional quality selections applied in our work. The profile of NGC 7099 is instead completely different, but the PMs used by Watkins et al. (2015a) were computed with a temporal baseline of only two years (see Bellini et al. 2014), and their quality is worse than that of our PMs.
The velocity-dispersion profile of NGC 6441 was also studied by Häberle et al. (2021). The PMs computed by Häberle et al. (2021) are obtained from a combination of data taken with ACS/HRC@HST and NACO@VLT detectors, which are better suited for probing the centermost arcsec of this very-crowded cluster. Our PM profile seems to suggest a more moderate increase of the velocity dispersion toward the center, although the values in the region in common within 10 arcsec from the center are in agreement at the 1σ level. The central velocity dispersion inferred in our work with the polynomial fit is σ r=0 µ = (0.285 ± 0.012) mas yr −1 , while their σ µ measured at r ∼ 0.76 arcsec is (0.316 ± 0.034) mas yr −1 . It is hard to clarify the nature of the discrepancy between these profiles given the different resolutions of the instruments used to infer the PMs in the centermost arcsec of the cluster.
Finally, we compared our PM velocity dispersions with those based on LOS velocities in the Baumgardt database, which are taken from various sources in the literature. There is an overall agreement between the PM-and LOS-velocity-based profiles, although there are discrepancies in some cases (NGC 5053,NGC 5272,NGC 5286,NGC 5466,NGC 6093,NGC 6101,NGC 6144,NGC 6205,NGC 6584,NGC 6681,NGC 6809,NGC 6981 and NGC 7099) where our PMbased σ µ are higher than the σ LOS , either in general or only in the centermost bins. The origin of the differences between PM-and LOS-velocity-based profiles might be related to systematics in either data sets, errors in the cluster distance (see Sect. 5.1), or instead be a proxy of a peculiar kinematic state of the cluster. However, a detailed comparison between these σ µ is outside the scope of this paper.

GENERAL KINEMATIC PROPERTIES
The collection of kinematic profiles in Figures A2-A5 allows us to analyze some of the general properties of stellar clusters, similar to what has already been done in the literature (e.g., Watkins et al. 2015a), but with a larger sample. We cross-correlated the kinematic pieces of information derived by means of HST and Gaia data with the structural properties of the GCs.
A quantity of particular interest in the kinematic characterization of star clusters is the anisotropy of the velocity distribution. Simulations following the evolution of star clusters during the initial violent-relaxation phase and including the effects of the tidal field of the host galaxy have shown that these systems emerge from this early evolutionary phase with an isotropic velocity distribution in the core, a radially anisotropic distribution in the intermediate regions, and an isotropic or slightly tangentially anisotropic distribution in the outermost regions (see Vesperini et al. 2014). This specific configuration depends on the initial conditions at the formation of the cluster, on the surrounding environment, and it changes during the subsequent long-term evolution; in particular during a cluster's long-term evolution, the effects of two-body relaxation and mass loss lead to a gradual decrease of the radial anisotropy imprinted during the early dynamical phases (Tiongco et al. 2016b, see that paper also for the possible development of radial anisotropy in tidally underfilling clusters with an initial isotropic velocity distribution). Figure 6 shows the ratio between the tangential and radial components of the velocity dispersions measured at the half-light radius as a function of the dynamical ages of the clusters as measured by the ratio of their Figure 6. Anisotropy at the half-light radius as a function of the ratio of the cluster age to the half-mass relaxation time. An isotropic system is characterized by σtan/σ rad = 1. Red points refer to core-collapsed GCs. All other systems are shown as black dots. The three solid lines show the time evolution of σtan/σ rad from three Monte Carlo simulations with initial W0 and filling factors (defined as the ratio of the half-mass to tidal radius) equal to (5, 0.09; green line), (5, 0.18; blue line), (7, 0.06; orange line); the three systems reach core collapse, respectively, at Age/t rh 10.9, 24.7, 6.8.
physical age 13 to their half-mass relaxation time (t rh ; from the Baumgardt catalog 14 ). The values of σ tan and σ rad at the half-light radius were derived by fitting a 4-th order polynomial function to the corresponding profile as described in Sect. 4. We excluded from the plot all GCs for which the HST data do not reach the half-light radius (so the anisotropy ratio would need to be extrapolated). Red points represent core-collapsed GCs, while all other systems are presented as black points. Core-collapsed GCs (regardless of whether they are considered as 'possible' or 'post' core-collapsed GCs) were labeled as such according to Trager et al. (1995)   From top to bottom, the plots refer to clusters with Age/t rh ≥ 10, 7 ≤ Age/t rh < 10 and Age/t rh < 7. Left panels refer to non-core-collapsed clusters, while those on the right shows the results for core-collapsed clusters.
Harris catalog 15 . Our analysis shows that dynamically older clusters (Age/t rh 10) tend to be isotropic even at the half-light radius, while dynamically-young systems are characterized by a radially-anisotropic velocity distribution at the half-light radius. A transition between the two regimes happens at Age/t rh between 7 and 10. These findings are consistent with the theoretical expectations discussed above (see, e.g., Vesperini et al. 2014;Tiongco et al. 2016b;Bianchini et al. 2017) and previous observational works (Watkins et al. 2015a).
To further illustrate the theoretical expectations concerning the evolution of the radial anisotropy, in Fig. 6 we show the time evolution of σ tan /σ rad (calculated at 15 There are some discrepancies between the list of core-collapsed GCs in Trager et al. (1995) and those in the Harris catalog, specifically NGC 6717 and NGC 6723. In the following, we consider the classification of Trager et al. (1995), for which NGC 6717 is a possible core-collapsed GC and NGC 6723 is not.
the projected half-light radius) from a few Monte Carlo simulations run with the MOCCA code (Giersz et al. 2013). The simulations follow the dynamical evolution of a few simple stellar systems composed of 500k stars with masses following a Kroupa (2001) initial mass function between 0.1 and 0.8 M , and spatially distributed according to the density profiles of King (1966) models with values of the central dimensionless potential equal to W 0 = 5 and W 0 = 7 (corresponding, respectively, to c 1.03 and c 1.53). The systems are characterized by an initial anisotropic velocity distribution following the Osipkov-Merrit profile (see, e.g., Binney & Tremaine 2008), β = 1 − σ 2 tan /(2σ 2 rad ) = 1/(1 + r 2 a /r 2 ) with r a equal to the half-mass radius. As shown in this plot, the initial anisotropy of the clusters at the half-light radius gradually decreases during the cluster's long-term evolution. For the tidally filling system, the enhanced rate of star loss leads to a more rapid isotropization of the velocity dispersion. At the time when the system reaches core collapse, the cluster's radial anisotropy slightly increases, then continues its gradual decrease toward isotropy (see also Tiongco et al. 2016b, for a study of the evolution of anisotropy for systems with a variety of different initial conditions). The differences between the anisotropy of systems that have similar dynamical ages but are in the pre-or post-core-collapsed phase is small, and within the uncertainty of the observed values.
In Fig. 7, to further investigate the kinematic anisotropy in stellar clusters, we divided the sample into three groups, i.e., clusters with a Age/t rh ≥ 10, between 7 and 10 (the Age/t rh transition region found in Fig. 6), and lower than 7. In each group, we also separated core-collapsed clusters (right panels) from all other systems (left panels). We collected all anisotropy measurements shown in Figs. A2-A5, and plotted them as a function of distance from the center of the cluster in units of r h . Gray points are the individual measurements, while black dots are the moving averages of those points. Clusters with Age/t rh ≥ 10 are isotropic at all distances within our FoV, regardless of their corecollapsed status, which is what we expected. Clusters with Age/t rh between 7 and 10 are again shown to be isotropic at all distances, although the non-corecollapsed sample hints at a marginal radial anisotropy at r r h . Finally, dynamically-young clusters clearly present the expected radial anisotropy at large radii.
Thus, clusters that underwent a core collapse seem to have similar velocity fields as those of the other GCs with similar dynamical ages; this appears to be consistent with what suggested by the results of the simulations presented in Fig. 6, which show that core collapse has only a relatively small effect on the radial anisotropy measured at the half-light radius. Additional simulations and a larger observational sample of core-collapsed clusters are necessary to further explore this issue. In particular, it is worth noticing that the core-collapsed sample with Age/t rh < 7 is composed of only four clusters: NGC 6541, NGC 6752, NGC 7099 and NGC 7078. While the first three objects show an isotropic field even slightly farther than the half-light radius, the latter presents a strong radial anisotropy. This feature for NGC 7078 has also been noted by Bellini et al. (2014) and . Among these four GCs, NGC 7078 (i) is more massive, (ii) is further from the center of the Galaxy, and (iii) has had fewer interactions with the Galactic potential of the Bulge/Bar and Disk (see . The combination of these properties could have preserved some of the original radial anisotropy in the innermost regions (e.g., Vesperini et al. 2014). Core-collapsed GCs are not only more spatially concentrated than other GCs, but their velocity-dispersion radial profiles are also steeper. This has been shown by Watkins et al. (2015a) and Cohen et al. (2021). In Fig. 8, we provide an updated version of this finding, and a comparison with theoretical models. The bottom panel of Fig. 8 shows the ratio between the velocity dispersion at the core and half-light radii (σ rc /σ r h ), computed as described before at r c and r h from the Baumgardt catalog, as a function of the concentration index c. The values of c are taken from McLaughlin & van der Marel (2005) and are defined as the log(r t /r 0 ), where r t and r 0 are the tidal and King scale radii, respectively. We selected c as obtained from the fit of a King (1966) profile. If c is not provided for a GC in the work of McLaughlin & van der Marel (2005), we used the value in Harris (1996Harris ( , 2010. Red, filled circles mark core-collapsed GCs, while all other systems are plotted as black dots. Green points are clusters taken from the sample of Cohen et al. (2021). This panel highlights a different location in the plot for core-collapsed GCs and the other stellar systems. This Figure confirms the trend between the spatial concentration and the steepness of the radial profile of the velocity dispersion.
In order to further explore this trend and carry out a consistent comparison with theoretical models, we show in the top panel of Fig. 8 Fig. 8), which, according to those authors, are in an advanced dynamical state, being close to (or having recently undergone) core collapse. Although this issue requires further investigation, it may represent a kinematic manifestation of the deviation of clusters from the dynamical properties of King models as they approach core-collapse.
There are a few exceptions in the overall picture described in the top panel of Fig. 8. The core-collapsed GC with c < 2.0 and σ rc /σ r h < 1.2 is NGC 362. NGC 362 is a post-core-collapsed GC (e.g., Dalessandro et al. 2013;Libralato et al. 2018a), and its peculiar position could be explained by the structural and dynamical evolution of post-core-collapsed GCs being driven by gravothermal oscillations (e.g., Makino 1996).
Three clusters shown with black dots seem to have characteristics similar to core-collapsed GCs. These GCs are NGC 5272, NGC 6652 and NGC 6715. McLaughlin & van der Marel (2005) pointed out that NGC 5272 shows deviations from the classical King isotropic model, and can be better fit with other models (see also Da Costa & Freeman 1976;Gunn & Griffin 1979). NGC 6652 is not considered a core-collapsed GC, but it is very concentrated and has a very steep surface-brightness profile, typical of core-collapsed systems (Noyola & Gebhardt 2006). NGC 6715 is a cluster 16 https://github.com/mgieles/limepy at the center of the Sagittarius Dwarf, an environment that could explain its peculiar core-collapsed-like concentration.

Kinematic distances
The Baumgardt's GC database also contains LOS velocity measurements. We made use of them to estimate distances of GCs by using the simple relation between the velocity dispersion along the LOS and in the plane of the sky: σ LOS = 4.7404 d σ µ ; where d is the distance in kpc, σ LOS and σ µ are the velocity dispersions along the LOS (in km s −1 ) and in the plane of the sky (from PMs in mas yr −1 ), respectively, and 4.7404 (km yr kpc −1 mas −1 s −1 ) is the conversion factor. The only assumption here is that the GCs are isotropic. We computed the distance d by comparing the values of the velocity dispersion at the same distance, i.e., at the center of the cluster, obtained for both LOS and PM measurements by fitting a polynomial function to the corresponding velocity-dispersion radial profiles as described in Sect. 4. Our results are summarized in Table 2. We considered only clusters with enough LOS velocities to solve the polynomial fit. Figure 9 shows comparisons between our distances and those from Baumgardt & Vasiliev (2021, panel a;using Gaia parallaxes, the comparison between PM and LOS velocities, and/or star counts), the Harris catalog (panel b; collected from various sources in the literature), the work of Watkins et al. (2015b, panel c; obtained with an approach similar to that used in our paper) for a sample of 14 GCs, and the estimates from Recio-Blanco et al. (2005, panel d; using the luminosity level of the zero-age horizontal branch 17 ). The median differences between d in our work and those in the literature shown in Fig. 9 are (−0.01 ± 0.16) , (−0.05 ± 0.18) kpc (Harris catalog), (−0.19 ± 0.17) kpc (Watkins et al. 2015b), and (−0.23 ± 0.20) kpc (Recio-Blanco et al. 2005), respectively. All distance estimates are in agreement with these literature values at the ∼1σ level. This is further proof of the goodness of our PM measurements. At large distances (>10 kpc), the differences between our distances and those from the literature increase, and so does the scatter of the points in Fig. 9. The lower velocity dispersions of some of these clusters, as well as the larger uncertainties both in 17 Distance moduli in HST Wide-Field Planetary Camera 2 (WFPC2) F555W filter of Recio-Blanco et al. (2005) were converted to distances in kpc after correcting for extinction using the extinction coefficient A F555W provided in Holtzman et al. (1995) and the E(B − V ) reddening in the Harris catalog. the PMs, LOS velocities and parallaxes, are likely the reasons for these discrepancies.

POSSIBLE APPLICATIONS
We choose the GC NGC 5904 to showcase some scientific applications enabled by our PMs 18 . Note that not all PM catalogs provide the same PM precision and overall quality, and some of the examples described below cannot be applied.

Internal kinematics of mPOPs
Kinematic differences between mPOPs can be a proxy of the different initial formation and evolution of first-(1G) and second-generation (2G) stars (Bekki 2010;Mastrobuono-Battisti & Perets 2013, 2016Hénault-Brunet et al. 2015;Tiongco et al. 2016bTiongco et al. , 2019Vesperini et al. 2013Vesperini et al. , 2021Calura et al. 2019;Sollima 2021). Recently, numerous observational efforts have investigated the kinematic properties of 1G and 2G stars to help us understand the mPOP phenomenon. We now know that 1G and 2G stars share similar kinematic features 18 NGC 5904 shows one of the cleanest and best-defined rotation curves obtained for a GC so far , thus suggesting that no significant residual rotation affects its kinematics in the plane of the sky. in dynamically-old GCs, or at least in the centermost regions of GCs where the two-body relaxation time is short, because two-body processes have already erased any kinematic differences between mPOPs (Anderson & van der Marel 2010;Libralato et al. 2018aLibralato et al. , 2019. The outermost regions of GCs are instead less relaxed, and some fingerprints of the initial kinematic differences between mPOPs can still be detected (Richer et al. 2013;Bellini et al. 2015Bellini et al. , 2018Milone et al. 2018;Dalessandro et al. 2018aDalessandro et al. ,b, 2019Cordoni et al. 2020a,b). Evidence of the presence of mPOPs in NGC 5904 has been shown both spectroscopically (Ivans et al. 2001;Carretta et al. 2009b;Gratton et al. 2013) and photometrically (Milone et al. 2017;Lee 2017Lee , 2021. To identify the mPOPs in our field, we made use of the pseudo color-color diagram ("chromosome map") computed by Milone et al. (2017) for this stellar cluster. We cross-correlated our PM catalog of NGC 5904 with that of Milone et al. (2017) and used their chromosome map to select 1G and 2G stars along the RGB of this cluster (right panels of Fig. 10). We identified three groups of stars: a 1G population (hereafter POPa in red) and two 2G groups (POPb and POPc in azure and green, respectively). The 1G-2G tagging was obtained in a similar way to what is shown in Fig. 4 of Milone et al. (2017). The two 2G sub-populations were arbitrarily identified using the Hess diagram in Fig. 10.  We then measured 19 the velocity dispersions of each population in various radial bins of at least 50 stars each. The velocity-dispersion (bottom-left panel) and anisotropy (top-left panel) radial profiles show that 1G and 2G stars have similar kinematic temperature and are isotropic within our FoV (there is only a marginal hint of a radially and tangentially anisotropic POPa and POPc, respectively). This is expected given that our field covers out to about the half-light radius of the cluster, a region where two-body encounters have likely removed any initial kinematic differences between mPOPs. Our findings are in agreement with those obtained by Cordoni et al. (2020a) with the Gaia DR2 PMs.
This result is in contrast with the finding of Lee (2021), which measured the analog of the POPc stars as more spatially concentrated than the other two populations even within the half-light radius. However, analogously to the internal motions, differences in the spatial segregation of mPOPs can be preserved in regions where the relaxation time is long enough to preserve them. Thus, the complete mixing and similar kinematic features are likely expected for the mPOPs in the core of NGC 5904.
To shed light on this disagreement, we computed the spatial distribution of the three populations in our FoV as follows.
First, we divided our sample 20 of RGB stars into three equally-populated bins and computed a kernel-density distribution of the ∆ * C F275W,F336W,F438W color for the stars in each bin. The ∆ C F275W,F336W,F438W color was obtained Milone et al. (2017) after rectifying the RGB sequences in the m F814W versus C F275W,F336W,F438W CMD. However, the sequences have been rectified using all stars in the FoV and, in different radial bins, they can show some deviations from being exactly vertical. Thus, we finetuned the ∆ C F275W,F336W,F438W color in each radial bin by using two fiducial lines (one each on the red and blue sides of the RGB sequence, respectively) drawn by hand and then computing a new color ∆ * 19 In addition to the criteria described in Sect. 4, stars that were analyzed in this section passed the photometric criteria described in Milone et al. (2017). 20 We removed the constraint on the PM error to increase the statistics at our disposal.
The kernel-density distribution 21 of the ∆ * C F275W,F336W,F438W color in each bin was fitted with a triple-Gaussian function (central panels of Fig. 11), and the fraction of stars in each population was estimated following the same approach of Bellini et al. (2013). To ensure the density distribution was robust against the small statistics along the RGB, we obtained the fraction of each mPOP by bootstrapping with replacements the sample of stars 1 000 times. The final values for the fractions of stars and their errors were determined as the median and the 68.27-th percentile of the distribution about the median, respectively. The ratios of the three mPOPs are shown in the right panel of Fig. 11. At odds with the finding of Lee (2021), we can see that the fractions of each mPOP do not vary as a function of distance from the center of the cluster within our FoV, as expected.

Energy equipartition
The energy-equipartition state of a GC is one of the most challenging measurements to obtain because it requires precise PMs along a (relatively) wide range of masses, i.e., for faint stars. Nowadays, this is one of the few applications that only HST can allow us to investigate in detail.
In collisional systems in a certain state of energy equipartition, there is a relation between the stellar mass m and the velocity dispersion σ µ : σ µ ∝ m −η , where η is the level of energy equipartition of the system. Theoretical works (Trenti & van der Marel 2013;Bianchini et al. 2016;Webb & Vesperini 2017) have shown that a complete state of energy equipartition (η = 0.5) is never reached because of the so-called Spitzer instability. Recent observational works on this topic have verified the goodness of this prediction (Anderson & van der Marel 2010;Bellini et al. 2018;Libralato et al. 2018aLibralato et al. , 2019. We made use of the exquisite PMs of NGC 5904 to measure its state of energy equipartition using both the parameter η defined above and the equipartition mass parameter m eq introduced by Bianchini et al. (2016) 22 . As a reference, a high degree of energy equipartition is characterized by a large value of η and a small value of m eq . We started by inferring the mass of stars along the MS below the MS turnoff by means of the updated isochrones of the "A Bag of Stellar Tracks and Isochrones" (BaSTI; Hidalgo et al. 2018). The parameters for NGC 5904 were chosen from the recent work of Gontcharov et al. (2019): a solar-scaled, 12.15 Gyr old isochrone for [Fe/H] = −1.33, Y = 0.2478, accounting for atomic diffusion, with a distance of 7.4 kpc. The fit is shown in Fig. 12. Although not perfect for the faintest portion of the MS due to the still-existing shortcomings in this mass regime of the color-effective temperature relationship, the fit is good enough to assign a mass to each star.
We then computed the σ µ in 35 bins of 1077 stars each along the MS. The bottom-left panel of Fig. 12 shows the result obtained by fitting the mass-dependent exponential relation of Bianchini et al. (2016) with a maximumlikelihood approach. We find a global m eq = 1.18 ± 0.07 M . Bianchini et al. (2016) provided a relation between m eq and the ratio between the cluster age and its core relaxation time. Using the age of Gontcharov et al. (2019) and a core relaxation time of 0.19 Gyr (Harris 1996(Harris , 2010 edition), we expect a value of m eq = 1.67 +0.32 −0.28 M , which is in agreement with our estimate at the 2σ level.
The middle-left panel presents the linear fit in the same log-log plane. We found a value of η = 0.25 ± 0.01. As for the value of m eq , this finding is consistent with the theoretical predictions of a partial state of energy equipartition even in an advanced stage of the cluster evolution.
As a reference, the median 1-D PM error as a function of stellar mass is plotted in the top-left panel of Fig. 12. The three lines correspond to the median PM trends of: all objects in the catalog (red line), sources that passed our photometric quality selections (see Sect. 4; yellow line), and stars that survived both the photometric and astrometric cuts (always described in Sect. 4; green line).
Our data set allowed us to push this investigation even further. We divided our FoV in four equally-populated 22 According to Bianchini et al. (2016), the relation between the velocity dispersion σµ and the stellar mass is as follows: radial bins, and measured the level of energy equipartition in each bin using the velocity dispersions measured in ten equally-populated magnitude bins. The result is shown in the bottom panels of Fig. 13 (left for m eq and right for η). We can notice marginal evidence of the innermost regions of the cluster having a higher degree of energy equipartition (low m eq and high η) than the outskirts. The innermost point within r c is the most uncertain because the fit is obtained with a smaller mass baseline.
In the top panels, we reproduce a similar analysis, but this time measuring the level of energy equipartition using the radial and tangential components of the velocity dispersion separately. We find that the levels of energy equipartition from σ rad and σ tan are consistent with each other at all radii, with only marginal differences at the 1-2σ level. The larger difference between the level of energy equipartition in the two components is shown in the innermost bin. The discrepancy is mainly the result of a poor fit, especially in the case of the exponential fit. This is in agreement with the simulations of Pavlík & Vesperini (2021, 2022 who found the degree of energy equipartition in the two velocity components to be similar in the region within the half-light radius. Future studies will extend the investigation of the energy equipartition in the tangential and radial dispersions to the outer regions of GCs where, according to Pavlík & Vesperini (2021, 2022, the degrees of energy equipartition in these two velocity components may differ. A note on the PM errors for faint stars -Despite our careful data reduction and PM computation, the PM errors might be under/overestimated. For example, a higher velocity dispersion for faint stars could be caused by underestimated PM uncertainties. We repeated the analysis shown in Fig. 12 by using only stars brighter than (1) m F606W = 22.1 (∼ 0.5M ) and (2) m F606W = 23.2 (∼ 0.42M ). We find: • m F606W ≤ 22.1 (M 0.5M ): m eq = 1.60 ± 0.18 M , η = 0.21 ± 0.02 ; • m F606W ≤ 23.2 (M 0.42M ): m eq = 1.31 ± 0.10 M , η = 0.24 ± 0.02 .
Although the results are in agreement at the 1-2σ level, it seems that the value of m eq is larger (η is smaller) the brighter is the magnitude cut, i.e., when we exclude the faintest objects with large PM errors. However, the exponential fit is less robust when the mass interval is small (the associated errors increase). The variation of η is instead mild.
Furthermore, the selections applied can bias the analyzed sample and the inferred kinematics, as extensively discussed in Bellini et al. (2014). In the case of   NGC 5904, we find that the selection on the PM reduced χ 2 applies a more significant cut on the PM errors, as shown in the top panel of Fig. 12 by the comparison of the yellow and green lines. As a test, we repeated the same analysis as before upon rescaling the PM errors of the stars used in the measurement of the level of energy equipartition by the ratio between the median PM errors obtained with sources that passed our photometric quality selections (yellow line in the top panel of Fig. 12), and with stars that survived both the photometric and astrometric cuts (green line). We find a global level of energy equipartition of m eq = 1.50 ± 0.11 M and η = 0.20 ± 0.01. Since we increased the PM errors, the intrinsic velocity dispersion of faint stars decreased, but again all estimates are in agreement within 2-3σ.
It is hard to say if these tests are a proxy of underestimated PM errors for faint stars,if they are biased by the less-constrained fit, or if what we see is just due to the intrinsic kinematics of our tracers. Nevertheless, these examples highlight how important it is to understand the data set used. Once again, we advise users to carefully check the PMs, and test their quality selections on a cluster-by-cluster basis, especially when the PM errors are of the order of the intrinsic velocity dispersion of the stars.

Kinematic mass determination
Knowledge of the relation between mass and velocity dispersion allows us to measure the kinematic mass of stars, similar to what has been done by Baldwin et al. (2016) and Libralato et al. (2018aLibralato et al. ( , 2019. The unknown mass of a group of stars X (M X ) with velocity dispersion σ X can be obtained as follows: where M ref and σ ref are the known mass and velocity dispersion of a reference population. For this analysis, we choose to measure the mass of two groups of stars: blue stragglers (BSs) and white dwarfs (WDs). Our reference population is the stars brighter than MS turnoff along the SGB and RGB (the same stars analyzed in Fig. A3), for which we assume M ref to be equal to the mass of stars at the MS turnoff. Using the isochrone in Sect. 6.2, we define M ref = 0.78 M . We selected BSs and WDs in various optical and UV CMDs (right panel of Fig. 14). We split the BSs into two groups (of 22 and 21 stars, respectively) and WDs into one group (of 14 stars), and measured their velocity dispersions σ µ . The result is shown in the bottom-left panel of Fig. 14. Blue points depict the velocity disper-sions of the BSs, while the green dot refers to the WD kinematics. Black points are the reference population. The effect of the energy equipartition is shown in the plot: BSs, which are more massive than RGB, SGB and MS stars, seem to be slightly kinematically colder (lower σ µ ) than the reference population. WDs have a mass of about 0.5 M (e.g., Bedin et al. 2019, and references therein). Being less massive than the reference population, their σ µ is instead marginally higher. Finally, the top-left panel shows that BSs and WDs are kinematically isotropic as the other stars at the ∼1σ level.
We previously fit the σ µ of the reference population with a polynomial function (see Sect. 4). To obtain the mass of the BSs, we repeated the same computation but we also solved for α BS by rescaling the polynomial of the reference population to fit the BS σ µ at the same time. By assuming the global value of η = 0.25 ± 0.01, we obtain M BS = 0.84±0.28 M . However, we notice that the BSs are preferentially located in regions where the local level of energy equipartition is higher than the global level. If we use the average value of η in these regions (0.34 ± 0.06), the mass becomes M BS = 0.82 ± 0.20 M . The difference between the two values is small since the value of α BS is close to unity. Both estimates are in agreement with the BS mass of NGC 5904 obtained by Baldwin et al. (2016, 0.82 +0.29 −0.18 M ). The mass of the BSs in this cluster is slightly lower than the average BS mass (1.0-1.6 M ; e.g., Ferraro et al. 2018). However, it is worth noticing that our sample is mainly composed of BSs close to the MS turnoff, which are less massive than those on the bright-end of the BS sequence (e.g., Raso et al. 2019).
For the WDs, a robust analysis cannot be obtained because we have only one point at our disposal for deriving α WD . Nevertheless, we repeated the same analysis to qualitatively assess the mass of the WDs. If we use the global value of η, we obtain M WD = 0.39 ± 0.25 M , which is lower than the average mass of the WDs in GCs. However, if we assume the local level of energy equipartition at the average WD distance from the center of the cluster, we find M WD = 0.48 ± 0.21 M .

CONCLUSIONS
We computed PMs for 57 stellar clusters studied in the GO-13297 program. The astro-photometric catalogs taht we have made publicly available represent the most complete, homogeneous collection of PMs of stars in the cores of stellar clusters to date, and more than double the number of clusters for which highprecision, HST -based PMs are available (Bellini et al. 2014). Furthermore, the astrometric information that we are releasing is complementary to that provided by the current (and future) Gaia data releases. At the dawn of a new era in astronomy with the first light of the James Webb Space Telescope (JWST ), the legacy that these PM catalogs offer is further enhanced, since they can serve as an important astrometric benchmark for JWST -based data reduction and tools.
We described the data reduction and, in great detail, the quality selections needed to select reliable objects for any kinematic analysis. We stress again that the data used for each cluster are different, thus any correction or selection should be tailored on a cluster-by-cluster basis. This is particularly important for stars with PM errors of the same order as the amplitude of the kinematic features that one wants to measure, for example, for stars along the MS of GCs.
We made use of our catalogs to study the general kinematic properties of the bright, massive stars in our clusters. We provided additional evidence supporting early findings that dynamically-young systems have a radially anisotropic velocity distribution at the half-light radius, while in dynamically older clusters the velocity distribution is isotropic at the same distance from the center of the cluster. This trend is consistent with the theoretical results of simulations showing that initially radially anisotropic clusters evolve toward an isotropic velocity distribution during their long-term evolution. Interestingly, core-collapsed clusters show similar properties to the non-core-collapsed systems, although a larger sample of core-collapsed GCs will be necessary to confirm the similarities with non-core-collapsed clusters (in particular for the group with longer relaxation times).
Finally, we showcased our PM catalogs using the GC NGC 5904. We separated the mPOPs along the RGB of the cluster and showed that, within our FoV, 1G and 2G stars have the same kinematics, are kinematically isotropic and have the same flat radial distributions. A detailed analysis of the kinematics of mPOPs will be the subject of another paper in this series. We investigated in detail the level of energy equipartition of NGC 5904. This cluster is in an advanced stage of its dynamical evolution, yet it has reached only a partial state of energy equipartition, as predicted by theoretical simulations. Knowledge of the level of energy equipartition also allowed us to measure the kinematic masses the BSs and WDs, finding a good agreement with the typical masses of these objects obtained with different methods in the literature.
We cross-correlated each of our PM catalogs with the Gaia EDR3 catalog and computed the PM zero-points to transform our relative PMs onto an absolute system.
We considered only cluster members with wellmeasured HST PMs (see Sect. 4) and whose Gaia EDR3 PMs have an RUWE better than 1.25, an astrometric excess noise less than 0.4, a number of bad along-scan observations less than 1.5% of the total number of along-scan observations, a PM error in each coordinate better than 0.1 mas yr −1 , and G > 13. If fewer than 25 objects were found in common with our HST sample, we relaxed these parameters to increase the statistics. The PM zero-point in each coordinate is defined as the 3.5σ-clipped average value of the difference between the HST and Gaia PMs. We set the error equal to the error of the mean. We also added in quadrature to our uncertainties a systematic error for the Gaia EDR3 PMs of 0.026 mas yr −1 (obtained using Eq. 2 of , assuming an angular separation θ = 0 deg since we are analyzing the clusters' cores). We do not take into account the rotation of the clusters in the plane of the sky, which is included in the PMs from the Gaia catalog but not in those made with the HST data (Sect. 3), although the scatter due to this effect is small. All values are reported in Table A1.
This astrometric registration allows us to directly compare our PMs with those in the Gaia EDR3 catalog. We show the result in Fig. A1. The red lines are the plane bisectors, not a fit to the data. The tight alignment of the points to the plane bisectors shows that our PMs are consistent with those of the Gaia EDR3 catalog, and that our absolute registration is accurate. The large scatters along the y directions are due to the poor quality of the Gaia PMs in some clusters, likely because of crowding.

B. DESCRIPTION OF THE PUBLICLY AVAILABLE CATALOGS
We release the astro-photometric catalogs of the 57 stellar clusters of the GO-13297 program through the MAST archive 23 . Table A2 presents the columns of our PM catalogs. For each cluster, we also release the photometric catalogs obtained with the second-pass photometry discussed in Sect. 2 for each filter/camera/epoch the data of which were used to compute the PMs. These catalogs will allow users to reproduce the quality selections we applied in Sect. 4. The description of a typical photometric catalog is provided in Table A3. When using our catalogs, users might notice peculiar features in some CMDs. For example, there are some filter combinations that show a clear, yet nonphysical, split in the SGB and RGB. The sources in these anomalous branches are close to the saturation and present 23 DOI: 10.17909/jpfd-2m08.
See also: https://archive.stsci.edu/hlsp/hacks. large QFIT and RADXS values, so they can be easily removed from the analysis. Our photometric catalogs can contain saturated stars for a given epoch/camera/filter if these sources are not saturated in at least two other epoch/camera/filter combinations so as to enable a PM measurement. As discussed in Bellini et al. (2017a), KS2 does not deal with saturated pixels, and photometry for saturated stars is instead provided by the first-pass stage of data reduction. Even though the photometric systems of saturated and unsaturated stars were designed to be the same, we have sometimes noticed differences between the regions in the CMDs dominated by objects belonging to either of the groups. These differences can be small zero-points (thus creating a discrete discontinuity in the CMD) or more complex behaviors (like a spread in the CMD or a different RGB slope with Figure A1. Comparison between the HST (including the PM zero-points in Table A1) and Gaia PMs. The red lines are the plane bisectors, not a fit to the data. respect to that of the unsaturated objects). Caution is again advised when dealing with saturated sources.
The photometry of the same camera/filter at different epochs is registered on to the same VEGA-mag system. Small zero-point variations can still be present, but they are expected to be small (∼0.01-0.02 mag, i.e., of the order of the uncertainty in the VEGA-mag calibration).
The astro-photometric catalogs of NGC 362 were made public by Libralato et al. (2018a). We include these same catalogs in our online repository, and refer readers to the related paper for their description. For NGC 6352, we provide the astro-photometric catalogs used in Libralato et al. (2019).

C. KINEMATIC PROFILES AND TABLE
In this Section, we present velocity-dispersion and anisotropy radial profiles. We refer to Sect. 4 for the detailed description of Figs. A2, A3, A4, and A5. Table A4 provides the values of σ µ at r = 0 arcsec, r c and r h . The individual measurements are available on our website. Note. (i) The ID number of a source is the same in all catalogs for the same cluster. (ii) (X,Y) positions are defined at a specific epoch about halfway between the first and last epochs used in the PM computation. This reference epoch is provided in the header of the file. (iii) "Corr Flag" is set to 1 if the PM was a-posteriori corrected for high-frequency systematics, and to 0 if otherwise. Table A3. Description of a photometric catalog for one filter/camera/epoch. Bellini, A., Piotto, G., Milone, A. P., et al. 2013, ApJ, 765, 32