Hunting for intermediate-mass black holes in globular clusters: an astrometric study of NGC 6441

We present an astrometric study of the proper motions (PMs) in the core of the globular cluster NGC 6441. The core of this cluster has a high density and observations with current instrumentation are very challenging. We combine ground-based, high-angular-resolution NACO@VLT images with Hubble Space Telescope ACS/HRC data and measure PMs with a temporal baseline of 15 yr for about 1400 stars in the centermost 15 arcseconds of the cluster. We reach a PM precision of $\sim$30 $\mu$as yr$^{-1}$ for bright, well-measured stars. Our results for the velocity dispersion are in good agreement with other studies and extend already-existing analyses of the stellar kinematics of NGC 6441 to its centermost region never probed before. In the innermost arcsecond of the cluster, we measure a velocity dispersion of (19.1 $\pm$ 2.0) km s$^{-1}$ for evolved stars. Because of its high mass, NGC 6441 is a promising candidate for harbouring an intermediate-mass black hole (IMBH). We combine our measurements with additional data from the literature and compute dynamical models of the cluster. We find an upper limit of $M_{\rm IMBH}<1.32 \times 10^4\,\textrm{M}_\odot$ but we can neither confirm nor rule out its presence. We also refine the dynamical distance of the cluster to $12.74^{+0.16}_{-0.15}$ kpc. Although the hunt for an IMBH in NGC 6441 is not yet concluded, our results show how future observations with extremely-large telescopes will benefit from the long temporal baseline offered by existing high-angular-resolution data.


INTRODUCTION
Globular clusters are the oldest surviving stellar systems in the Galaxy. Because of their long dynamical timescales, they are witnesses of the early history of the Milky Way. The study of the internal dynamics of globular clusters is an active field of research and can reveal information about the formation and evolution of the clusters themselves, their interaction with the Galactic potential, but also the possible existence of intermediate-mass black holes (IMBHs) in their core.
IMBHs have been predicted in the centres of globular clusters (Portegies Zwart & McMillan 2002;Miller & Hamilton 2002), but a definitive detection is still lacking (see Greene et al. 2019 for a recent review). As the sphere of influence of a hypothetical IMBH is limited to the very centre of the cluster, dynamical studies of this central region are necessary to infer the presence of an IMBH.
There are two observational methods to study the kinematic signature of individual stars in a globular cluster: spectroscopic line-of-★ E-mail: haeberle@mpia.de sight (LOS) velocity measurements and astrometric measurements of proper motions (PMs). While in the early days of the field the line-of-sight velocities of only small numbers of stars were measured, in recent years significant progresses have been made both in spectroscopy and astrometry. The MUSE integral field spectrograph allowed LOS velocity measurements of up to 20,000 stars within the half-light radius of a large sample of globular clusters (Kamann et al. 2018).
On the astrometry side, the best tool to study the crowded cores of the clusters is the Hubble Space Telescope (HST) and there are catalogues with high precision PM measurements for up to a hundred thousand stars Libralato et al. 2018) in a single cluster. Additionally, the Gaia satellite (Gaia Collaboration et al. 2016) allows the study of stellar motions in the outskirts of globular clusters. However, its completeness and precision are limited in cluster centres, where the stellar density is high.
Observations of the crowded cores of some globular clusters still remain very challenging with current instrumentation. Crowding effects limit the usability of spectroscopic facilities, therefore, astrometry with high resolution imagers is the method of choice. One such example is the globular cluster NGC 6441: its core is extremely crowded and neither of the current HST imagers provides the necessary resolution for astrometric studies of the core. The Advanced Camera for Surveys Wide Field Channel (ACS/WFC) has a pixel scale of 50 mas pixel −1 and the Wide Field Camera 3 (WFC3) UVIS has a pixel scale of 40 mas pixel −1 .
The High Resolution Channel of the Advanced Camera for Surveys (ACS/HRC) had a smaller pixel scale of 25 mas pixel −1 and, while still affected by the high stellar density, could measure stellar position in the core of NGC 6441 with high precision. A first epoch of observations was taken in 2003, but the failure of this instrument in 2006 made it impossible to observe the core of NGC 6441 again.
Since there are no other suitable epochs of HST observations, we re-observed the core of the cluster using the near infrared instrument NACO at the Very Large Telescope (VLT) of the European Southern Observatory (ESO). The adaptive-optics (AO)-assisted observations close to the diffraction limit with an 8-m class telescope in the S band have a point spread function (PSF) full width half maximum (FWHM) of about 73 mas, similar to that of the HST at visible wavelengths, and with 13.2 mas pixel −1 the pixel scale is twice as small as the one of ACS/HRC. NGC 6441 is a Galactic bulge cluster. With a mass of 1.2 × 10 6 M (Baumgardt & Hilker 2018), it is one of the most massive clusters in the Galaxy. Furthermore, it is metal rich ([Fe/H] = −0.55, Harris 2010) and hosts at least two main stellar populations .
The dynamics of the cluster have been studied in several papers, both with PMs (Watkins et al. 2015a) and LOS velocities (Kamann et al. 2018), but due to its high density, accurate measurements of the velocity dispersion of the stars in the very centre of the cluster are still lacking. Our study probes the velocity dispersion to unprecedentedly small radii. While we cannot put strong constraints on the presence of an IMBH, we add more than 20 PM measurements in the innermost arcsecond. The combination of space-based and groundbased AO-assisted astrometric imaging with a long time baseline is a successful pilot study to showcase what will be possible with new instrumentation expected in the next decade, such as the MICADO imager at the Extremely Large Telescope (Davies et al. 2016).
This paper is divided in the following sections: first we describe our observations (Section 2), after which we present the data analysis in Section 3. The determination of the PMs is described in Section 4. We use our results to create dynamical models of the globular cluster in Section 5, and discuss and conclude our work in Section 6.

Epoch 1: 2003 HST ACS/HRC Observations
The HST instrument ACS/HRC was in operation from 2002 to 2006. It featured a 1024 × 1024 pixel CCD detector with a pixel scale of 28 × 25 mas 2 pixel −1 , giving it a field of view of ∼ 29 × 26 arcsec 2 .
The core of NGC 6441 was observed with the ACS/HRC in two broadband filters (F555W and F814W) during Program . The 36 exposures in the F555W filter have all the same exposure time of 240 s. The F814W band exposures comprise of 5 short and 12 long exposures of 40 s and 440 s, respectively.

The instrument
NACO (short for NAOS-CONICA) was a near-infrared, adaptiveoptics-assisted imager and spectrograph at the Very Large Telescope (VLT). A full description of the instrument and its performance can be found in Lenzen et al. (2003) and Rousset et al. (2003). The instrument was mounted on the Nasmyth B focus of UT4 from 2001 to 2013, and then moved to UT1 in 2014, where it continued its operations until its recent decommissioning in October 2019. The adaptive optics front end (NAOS) was equipped with a 185 actuator deformable mirror, a tip/tilt mirror and two wavefront sensors (operating in the visual and the IR range). Using bright natural guide stars, the best Strehl ratio obtainable in the S band was around 50% in typical observing conditions, similar to what we achieved in our observing run.
The camera CONICA was equipped with an Aladdin2 detector (1026 x 1024 pixels, InSb) that replaced the Aladdin3 detector that was in use from 2004 to 2013. The detector was affected by several artefacts that are described in detail later. CONICA offered cameras with three different pixel scales (S13: 13.22 mas pixel −1 , S27: 27.06 mas pixel −1 , S54: 54.3 mas pixel −1 ). For our observations of the core of NGC 6441, we made use of the S13 mode to obtain the highest possible resolution.

The dataset
The second epoch of observations of the core of NGC 6441 were obtained with the ESO Program ID 0101.D-0385 (PI: M. Libralato). We only made use of images taken during the August run because of their overall better quality. In these nights, there were 23 usable exposures of 200 s each, and 27 short observations of 30 s each. The observations were executed with a dither pattern that covered the field of the first-epoch HST observation almost completely (see Figure 1). The observing strategy was designed to solve for the geometric distortion of the detector using an auto-calibration approach (see, e.g., Libralato et al. 2014) and to achieve an astrometric precision high enough to allow for the kinematic analysis of NGC 6441.

DATA ANALYSIS
PM measurements require multiple epochs of precise stellar position measurements. In this section we describe how we extracted stellar positions from the raw exposures of our two datasets.

Epoch
Telescope of the HST PSFs from one exposure to the next even within the same telescope orbit.
Preliminary stellar positions and fluxes of bright sources were obtained through PSF fitting using the FORTRAN code hst1pass (Anderson in preparation, see Bellini et al. 2018c for details). Photometry of saturated stars includes all the relevant flux that has bled into adjacent pixels following the prescriptions given in Gilliland et al. (2010). Stellar positions were corrected for the effects of geometric distortion using the state-of-the-art solutions provided by Anderson & King (2004).
Next, we made use of Gaia DR2 positions (Gaia Collaboration et al. 2016 to define a reference-frame system with North up, East to the left, and with the same pixel scale of NACO of 13.2 mas pixel −1 . We transformed stellar positions of each singleexposure catalogue on to the reference frame by means of general, sixparameter linear transformations. Our best estimate of positions and fluxes for all possible sources in the ACS/HRC field is obtained using the FORTRAN code KS2 (Anderson in preparation, see Bellini et al. 2017a for details). KS2 starts from the image-tailored PSF models, the lists of bright stars and their transformations on to the reference frame, and goes through several waves of source finding, measuring and subtraction using all the exposures simultaneously. Our final firstepoch catalogue contains around 44 000 sources measured in both F555W and F814W filters. In addition to the photometry results, the output contains several quality parameters such as the radial excess value (RADXS). If it is positive, the profile of a single star contains excessive flux outside of the fit radius with respect to the PSF, if it is lower than 0, this flux is lower than expected.
To calibrate the photometry of the ACS/HRC data, we performed aperture photometry on the corresponding _drz images that are resampled and normalised to 1 s exposure time, corrected the results for finite aperture (using the corrections from Bohlin 2016) and then brought them onto the VegaMAG system using the zeropoints available at the STScI website 3 . Then we determined the zeropoint between the calibrated aperture photometry and our PSF photometry by taking the 3 -clipped median of the magnitude difference for bright isolated stars (we chose stars that had no brighter neighbours within a 18 pixel radius). A similar calibration process for HST photometry has been described e.g. in Bellini et al. (2017a).

Pre-reduction
We downloaded the raw NACO exposures and the corresponding calibration frames (dark frames, flat-fields, bad-pixel maps) from the ESO Science archive. The dark frames were not usable because the background showed patterns that varied over the course of the night. We therefore only divided the images by the flat-fields. Furthermore, we flagged all saturated pixels and added them to the bad-pixel map. The saturation threshold is set at 10 000 analogue-digital converter units (ADUs) to avoid non-linearity effects following the recommendations of the NACO User manual (Schmidtobreick et al. 2018). The bottom-left quadrant of NACO presents a high number of bad columns in a regular pattern (2 good, 2 bad, 3 good, 1 bad), so we chose not to use it in our analysis.

Background Model
Insufficient thermal shielding within the instrument 4 created a diffuse, non-static background pattern that cannot be corrected in the pre-reduction phase. The time-dependent nature of the background pattern and the high density of stars in the cluster centre made it challenging to correct this pattern. The correction we applied is the result of two iterations. In each iteration, we used the median of multiple images in which the background pattern did not change as our background model. To remove the influence of stars, we flagged pixels from an image if they were significantly brighter than in the other images. Furthermore, we flagged a circular area with = 45 NACO pixels around saturated stars to remove the influence of their bright and extended halos. Flagged pixels were excluded from the median calculation. In the first iteration this clipping was performed directly on the raw images. In the second iteration we improved the model by subtracting the modelled stellar fluxes from the raw images (based on the PSF photometry results from the next paragraph) before running a jackknife algorithm to find additional outliers. For some images it was beneficial to shift the model in the -direction before subtracting it. The different steps are visualised in Figure 2.

Determination of the PSF
An accurate PSF model is a crucial ingredient for high-precision astrometry. We created a 3 × 3 array of empirical PSFs to take into account the spatial variation of the PSF caused by the telescope optics and AO-related setup. Since the AO performance and seeing conditions change over time, we tailored these PSF arrays to each image. This method was introduced for ground-based instruments in Anderson et al. (2006) and has also been adapted to various near-infrared instruments (Libralato et al. 2014(Libralato et al. , 2015Kerber et al. 2019). We refer to these papers for a detailed description of the PSF modelling.
Here we provide only a brief overview of the method and discuss the major differences with respect to the original papers.
Because of the small field of view of NACO, we used a regular 3 × 3 grid of PSF models to map the spatial variability. Bilinear interpolation between the grid points was used to determine the PSF on each location of the detector. Due to the unusable lowerleft quadrant, the lower-left grid point contains no PSF. For every other gridpoint, the local PSF was determined from at least 20 wellmeasured, isolated bright stars in a nearby rectangular cell. The exact layout of these cells and the gridpoint locations are shown in Figure 3.
To determine the Strehl ratio of our observations, we divided the central value of our different PSF models by the central value of an Airy function with a radius of 69 mas (or 5.2 NACO pixels), the expected diffraction-limited PSF for a 8-m telescope observing at 2.2 µm wavelength. The Strehl ratios we obtained ranged between 0.2 and 0.49, with a median value of 0.37. This is compatible with the typical Strehl ratios stated in the NACO User manual. The median FWHM of our PSF was 73 mas. Figure 3 shows the 8 PSF models determined for an individual frame and how the local models differ from a spatially constant PSF model.

Fitting the stellar positions
After the PSF model has been determined for each image, we used it to fit the raw position (x,y) in pixels and the instrumental magnitude of (1,3) (2,3) (3,3) (1,2) (2,2) (3,2) (1,1) (2,1) (3,1) each star. To be able to select well measured stars, we also calculated the so called quality-of-fit (QFIT) value, which is a normalised sum of the fit residuals within the fitting radius (Anderson et al. 2006). The closer to 0 the QFIT, the better is the PSF fit.

Geometric distortion correction (GDC)
To achieve a sub-pixel astrometric precision, we have to correct the geometric distortion present in the NACO images, which reaches up to 2 NACO pixels in the corners of the detector. We redetermined the GDC using our HST catalogue as distortion-free reference, the process is described in detail in the Appendix A. We corrected the geometric distortion to a level of ≤ 0.03 NACO pixel (≈ 0.4 mas).
In addition to the GDC using an external reference, we also employed a so called autocalibration approach in which the distortion free reference is obtained from the NACO data itself. However, this lead to a worse correction. The failure of the autocalibration approach with our 2018 NACO dataset highlights the importance of the choice of the dither + rotation pattern, if such a calibration is required.

Creation of a NACO master frame and tests on the astrometric precision
Most stars in our dataset have been measured on multiple NACO exposures. To combine the distortion-corrected single-image catalogues, we matched them on the same reference system as defined for the first epoch (see 3.1). We transformed the stellar positions from each distortion-corrected single-image catalogue on to the reference system using six-parameter linear transformations. Stars that have been measured in at least three images were then added to our NACO master frame. Our best estimate for their position is the averaged position from the individual transformed positions.
The scatter of the single measurements, quantified by their root mean square (rms) deviation from the master frame position, is a measure of the astrometric precision we can reach for a given signal to noise ratio (S/N). In Figure 4, this positional rms is plotted as a function of the S magnitude (see subsection 3.2.7 for our flux calibration). As expected, we see a decrease in precision for faint stars, i.e., those with a lower S/N ratio. For stars brighter than S < 12.5, the trend flattens out at a level of approximately 0.03 NACO pixel (≈ 0.40 mas), while theoretical limits on centroiding precision (e.g. Lindegren 1978) predict: (with close to unity). This fundamental limit of 0.03 pixel is worse than the 0.01 pixel achieved for astrometry with HST instruments, but comparable to other ground-based IR studies (Libralato et al. 2014(Libralato et al. , 2015Kerber et al. 2019). It is most likely caused by systematic effects, such as residual uncorrected distortion, imperfections in our PSF models and/or systematic effects related to the thermal background pattern and its correction.

Photometric calibration and creation of colour-magnitude diagrams
We brought our NACO PSF photometry on to the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2006) photometric system by cross-matching it with the NGC 6441 catalogue published by Valenti et al. (2010) (in the following Val10). Due to the depth of the NACO images, most of the stars found in Val10 are saturated in our long exposures. Therefore, we used a two step-process. First, we created a catalogue based on the short NACO exposures whose photometry has been zeropointed on to our long-exposure master-frame. Then, we determined the magnitude difference between the 63 stars in common between our master frame and the Val10 catalogue. Our best estimate of the zeropoint is the sigma clipped median of the difference between the magnitudes in both catalogues. Since we only have one filter, we added this simple zeropoint without accounting for possible colour effects. After obtaining calibrated photometry for both the HST and the NACO dataset, we created ( F555W − F814W ) and ( F555W − S ) colour-magnitudes diagrams, which can be seen in Figure 5.

PROPER MOTIONS
To measure PMs one has to determine by how much the positions of individual stars have changed over time.
Compared to other studies, where multiple different epochs and fields of views are combined, the situation in our study is rather simple: we only combine two single epochs with a large time span in between. For this reason, we do not use a linear fit through multiple datapoints for each star, but we determine the position difference between two single epoch master frames. Due to the more reliably determined PSF and the better general astrometric precision, we decided to use only the long NACO exposures for this part of the analysis, at the cost of being unable to determine the PMs for around 20 stars, which were saturated in the long exposures.

Creating the single-epoch master frames
The goal of this step is to create new, improved master frames by using additional clipping procedures and local corrections.  steps were performed on both the HST ACS/HRC and the NACO data.
1. Selection of well-measured stars: we applied several selections based on quality criteria (QFIT value, RADXS value) to restrict our sample to well-measured stars. Those well-measured stars were then used to determine the parameters of the following transformations.
2. Global transformations: for each image, we use all stars flagged as "well-measured" to determine the optimal six-parameter linear transformations to transform the single-image catalogues on to the master frame. The residuals between the single-image catalogues and the master frame are used in the next steps for a local correction.
3. Local corrections: after the parameters of the global transformations have been determined, we used local corrections to remove residual local effects such as uncorrected distortion. For each star, we determined the clipped mean of the transformation residuals of the 50 closest neighbours in both the and the directions. These mean residuals are subtracted from the mean coordinates of the star. This procedure is called "boresight" correction (van der Marel & Anderson 2010).
4. Error-based clipping and determination of mean positions: the last steps gave us a list of multiple position measurements for each star, all transformed in the same reference frame. We defined as improved master-frame positions the mean positions of these locally-corrected measurements. Outliers were removed using a jackknife approach: for each star we excluded one measurement at a time and checked whether the excluded measurement deviates more than 10 standard deviations from the mean value of the remaining measurements. To determine the expected standard deviation, we employed the empirical error model based on the typical positional rms at a given magnitude (see Figure 4). After the clipping process, we calculated the mean of the positions based on the remaining measurements. As the single-epoch position error, we employed the standard error of the mean for both coordinates: Δ , master frame = , √

Combination of single-epoch master frames and determination of PMs
The procedures used to match the two single-epoch master frames are similar to the methods used for the master frame creation. The difference is that, instead of matching single-image catalogues from the same epoch, we now directly match the master frame positions of the two different epochs. Since we want to measure PMs relative to the bulk motion of the cluster, we use only bona-fide cluster members to determine the parameters of the linear transformations between the master frames.
1. Selection of well-measured cluster members: we use the same quality selection as described above. In addition, we included a selection of bona-fide cluster members based on their location on the ACS/HRC based ( F555W − F814W ) CMD (see Figure 5). Once PMs have been determined, we also restricted the sample of bonafide cluster members on the basis of their location on the vector-point diagram.
2. Global transformations: for each star, we use the positions of all well-measured stars within a radius of 1000 NACO pixels to determine the optimal six-parameter transformation between the two master frames. 5 3. Calculation of PM and PM error: the stellar PMs are obtained as the difference between the transformed positions of the two master frames, divided by the temporal baseline. The PM errors are calculated by quadratically adding the positional errors of the single epochs and then dividing this result by the temporal baseline.

A posteriori local corrections
By construction, the mean motion of cluster members should be at location (0,0) on the vector-point diagram regardless of stellar positions on the master frame. Local deviations of the mean motion are caused by small, uncorrected-for systematic effects, which we mitigated using an a-posteriori local correction as follows.
For each cluster star in our PM catalogue, we chose the 50 closest neighbouring cluster stars and calculated their mean motion. This value is subtracted from the measured PM of the star. This step leads to an additional statistical PM error of ∼0.028 mas yr −1 per coordinate. Therefore, we kept both the uncorrected and the corrected position residuals for the further analysis.

Resulting proper motion precision
The vector-point diagram and the PM uncertainties in both coordinates are shown in Figure 6. For bright stars with S < 13.5 (or F555W < 17.5) we reach uncertainties of around 0.03 mas yr −1 . Similar precisions are reached in pure HST based studies (see Bellini et al. 2014)  around 0.1 mas yr −1 for stars with G=17 mag even in less crowded fields (Lindegren et al. 2018).

Determination of the velocity dispersion
Only cluster stars (within 1.5 mas yr −1 from the origin of the vectorpoint diagram) that are well-measured in the HST and NACO data were used in the analysis of the internal kinematics. For the HST data, we defined as "well-measured" stars that have: (i) magnitude rms in both filters lower than 0.025 mag; (ii) QFIT parameter larger than 0.985; (iii) the absolute value of the shape parameter RADXS lower than 0.03. For the NACO data, we only required the 1D positional errors to be lower than 0.4 NACO pixel. Stars with a PM error in either direction greater than 0.15 mas yr −1 or larger than half the local velocity dispersion of the closest 50 cluster stars were excluded from the analysis. We tested different quality selections, and find negligible differences in the resulting velocity dispersions. The parameters described above provide a good compromise between including as many stars as possible (for better statistics in the kinematic analysis) and excluding poorly-measured objects. Out of the around 1400 stars with a PM measurement, ∼1200 were included in the kinematic analysis.
The velocity dispersions were obtained by subtracting in quadrature the PM errors from the observed scatter of the PMs (van der Marel & Anderson 2010). We analysed the combined ( ), radial ( R ), and tangential ( T ) velocity dispersions as a function of distance from the cluster's centre in: i) one radial bin with all stars within 1 arcsec from the cluster's centre (23 stars); ii) four equally-populated radial bins between 1 and 5 arcsec (107 stars each); iii) eight radial bins using all stars outside the centermost 5 arcsec (seven groups of 97 stars and one with 87).
The results are shown in Figure 7 and in Table B1 in the Appendix. Our data are very consistent with the PM results of Watkins et al. (2015a) and the MUSE LOS measurements of Kamann et al. (2018), which are also shown in Figure 7. However, we cannot reproduce the dip in velocity dispersion in the innermost MUSE data point, that may be caused by crowding effects (Alfaro-Cuello et al. 2020). In the innermost arcsecond we measure a combined velocity dispersion of (0.316 ± 0.034) mas yr −1 . If we employ the dynamical distance estimate from Section 5 ( = 12.74 kpc) this corresponds to (19.1 ± 2.0) km s −1 .

Search for fast-moving stars
Fast moving stars in the centre of NGC 6441 could be a possible signature of an IMBH (Drukier & Bailyn 2003). The lack of LOS velocities does not allow us to obtain a 3D kinematic picture of the core of this cluster. However, we can still investigate the presence and nature of fast-moving objects in the plane of the sky thanks to our high-precision PMs.
We initially searched for stars with a PM that indicates a velocity higher than the escape velocity of the cluster as a result of the interaction with a potential IMBH. We considered fast-moving stars those with total PM between 1.26 and 1.5 mas yr −1 . The lower limit corresponds to the escape velocity of esc =76 km s −1 reported by Baumgardt & Hilker (2018) for NGC 6441. We find 4 stars in this velocity range. Stars with a PM > 1.5 mas yr −1 are not considered in our analysis because they are either mismatches between the NACO Median PM error of bright stars (0.031 mas yr 1 ) Figure 6. The left panel shows the vector-point diagram with our measured PMs for NGC 6441. The sample is restricted to stars that have been used for the kinematic analysis. The two right panels show the magnitude dependence of the PM errors for the stars in the sample. The red line marks the median PM error for stars brighter than K S = 13.5 (∼ F555W = 17.5). One can see that for these stars, a PM precision of around 0.03 mas yr −1 per PM component is reached. and HST catalogues, or are located in the outskirts of our field of view, far away from the cluster centre.
For each high-PM star, we investigated if its PM vector suggests a past passage within 2.5 arcsec from the centre of the cluster (for a similar procedure see Libralato et al. 2021). This specific radius is the influence radius of an IMBH with a mass of BH = 1.3 × 10 4 M , i.e., the upper limit of the IMBH mass we computed in Sect. 5.3.2. We find only one high-PM star with a PM vector consistent with an ejection caused by the interaction with an IMBH in the core of the cluster.
We also performed a statistical analysis of the total PMs as done by  for the globular cluster Centauri. We divided our sample in 12 radial bins of 100 stars each, and determined various percentiles of the total PM distribution in each such bin (see Figure 8). If an IMBH is harboured in the core of NGC 6441, we would expect a higher number of high-velocity stars closer to the centre. However, we do not detect a significant difference between the distribution within the innermost 2 arcseconds and the distributions in the outer bins at > 8 arcsec.
We can conclude that, even though we do detect a single fastmoving star that is geometrically associated with the cluster centre, this is not enough to clearly indicate the presence of an IMBH. Deeper observations would be necessary to expand our brightness-limited sample and only the combination of radial velocity measurements with proper motions can give the full velocity vector.

Dynamical Models
Kinematics tell us how fast the stars are moving inside the cluster. Dynamical models offer us a way to determine the underlying physics that make the stars move as they do. We ran dynamical models of NGC 6441 to determine its mass profile, anisotropy profile, and dynamical distance. We also used our models to assess whether we could determine the existence and properties of a possible IMBH at the centre of the cluster.

Setup of models
We used a compilation of 4 kinematic datasets. We started with the NACO PM catalogue of well-measured stars that was used to determine the velocity dispersion profile of the cluster in Section 5.1. We augmented this with the catalogue of HST PMs from Watkins et al. (2015b), which probe further out in the cluster than the NACO dataset alone. To this, we added the LOS velocity dispersion profiles from Kamann et al. (2018) and Baumgardt & Hilker (2018). We used the surface brightness profile from Trager et al. (1995) as a proxy for the surface number density profile from which the kinematic samples were drawn.
For the dynamical models themselves, we used the spherical Jeans Anisotropic Multi-Gaussian-Expansion (JAM, MGE) models (Cappellari 2008(Cappellari , 2015. The methodology broadly follows that described in Section 5.2 and Appendix C of Hénault-Brunet et al. (2019), and we refer to that paper for details. Here we briefly summarise the analysis.
In the JAM models, the stellar and mass densities are treated as MGEs. Each component of the MGEs has a width and a weight. We assumed that the widths were the same for the stellar and mass MGEs but allowed the weights to vary independently so as to fit a variable mass-to-light ratio (M/L). By giving each component of the stellar MGE an anisotropy, we were also able to fit for a variable anisotropy profile. We used 5 Gaussian components to fit the cluster. We tried initially with 6 and found that this was too many, so reduced to 5 and found that this gave much improved and non-degenerate results. Altogether, this gave us 20 free parameters: 5 widths , 5 stellar density weights , 5 mass density weights , and 5 anisotropies . There were two further parameters -the distance , and the IMBH mass BH -taking the total number of free parameters to 22.
In the spherical JAM models, the mean velocities are everywhere zero, and the velocity distributions are characterised by the dispersions in the projected-radial, projected-tangential, and LOS directions. The likelihood of the data given the model was calculated differently for the PMs and the LOS velocities due to the different datasets available. We treated the PMs discretely, that is we did not bin the stars, but instead for each star calculated the likelihood of observing a star with the measured velocity and uncertainty given the mean velocity and velocity dispersion predicted by the model at the position of the star, assuming Gaussian velocity distributions and uncertainties. For the LOS velocities, we had only binned velocity dispersion profiles, so we calculated the likelihood of measuring a given velocity dispersion and uncertainty given the velocity dispersion predicted by the model at the position of the bin.
The Gaussian widths, Gaussian weights and the IMBH mass were all fitted in log-space, as they had the potential to span many orders of magnitude. We used priors to restrict the range of certain parameters for both physical and computational reasons, but otherwise used flat priors within the allowed ranges.
To avoid degeneracies between different components and to ensure that the models fit for 5 separate components, we insisted that log +1 − log > 0.2. Additionally, we limited the widths of the innermost and outermost MGE components such that 1 > 25 and 5 < max / √ 3, where 25 and max are the projected radial po-sitions of the 25-th star in the PM dataset and outermost point in surface brightness profile, respectively. The first condition ensured that there were at least 25 stars inside of 1 to constrain the inner cluster properties, and the second condition ensured that the outer slope of the density profiles was at least 3, which is required for a finite system. The light and mass weights were assumed to be positive. We further added a lower limit on the M/L of each component such that / > 0.1, as values lower than this would be unrealistic physically. is a modified anisotropy where and are the velocity dispersions in the radial and tangential directions in a spherical coordinate system. This modified anisotropy is defined to be both symmetric about isotropy ( = 0) and finite in extent ( = ±1). In practice, we only sample between ≈ ±0.88 6 to avoid extreme anisotropies that are computationally intensive but not seen in real clusters.
Finally, we restricted how much the density and anisotropy profiles could change from component to component by insisting that −7 < Δ , /Δ < 3 and −2 < Δ /Δ < 2, where Δ = +1 − and similarly for Δ and Δ and Δ . In practice, very few components hit the limits of these ranges.
We set a lower limit of the central black hole mass of 0.1 M ; the BH mass was fit in log space and very small masses are indistinguishable in the model so this parameter has the potential to go to −∞ if not limited. No other restrictions were set. The distance was fit in linear space, and was assumed to be positive with no further restrictions.

Modelling results
We explored the available parameter space using the affine-invariant Markov Chain Monte Carlo software (Foreman-Mackey et al. 2013). We used 100 walkers and ran for 7 500 steps to be sure that the chains had fully converged. The resulting fits to the cluster are shown in Figure 9. In each panel, the coloured data points show the data with their uncertainties (in cases where there are no points, data is either not available or not measurable). The black solid lines show the median of the fits for the profile, the darker shaded region spans the 15.9 to 84.1 percentiles (equivalent to the 1-sigma confidence region for a Gaussian distribution) and the light shaded region spans the 2.5 to 97.5 percentiles (equivalent to the 2-sigma confidence region for a Gaussian distribution).
The left panels show from top to bottom the projected radial dispersion profile R , the projected tangential dispersion profile T , and the LOS velocity dispersion profile z . In the top and middle panels, the blue points show the profiles calculated from NACO PMs and the orange points show the profiles calculated from the HST PMs, however we stress that the fits were not done to these profiles but to the individual measurements. The profiles are shown here for visualisation purposes. In the bottom panel, the green points show the MUSE dispersion profile from Kamann et al. (2018) and the purple points the dispersion profile from Baumgardt & Hilker (2018); for these we did fit directly to the binned dispersion profiles. In general, these fits are very good. The models struggle in the outer regions where there are only LOS constraints but not PMs, and show some broadening in the centre where there are very few stars. The only poorly-fit point is the central point in the MUSE LOS dispersion profile, which the models were not able to constrain. Kamann et al. (2018) also point out this central dip in their dispersion profile; they note that it could potentially be due to crowding but argue that this is likely not the case as the cluster does not have a steep central surface brightness profile. That the models are unable to reconcile the PM and LOS dispersion profile may suggest otherwise and that the central MUSE data point is more affected by crowding than previously believed. This effect due to crowding is also clearly seen in the MUSE velocity dispersion profile of the cluster M54 (Alfaro-Cuello et al. 2019, 2020 and can be overcome by higher spatial resolution data (Alfaro-Cuello et al. in prep.).
The panels in the central column show the anisotropies, from top to bottom they are tangential over radial, radial PM over LOS, and tangential PM over LOS. For the top panel, we are able to show data (NACO PMs in blue and HST PMs in orange as for the PM dispersion panels) as we have the same set of bins for both the radial PM and tangential PM datasets. For the middle and lower panels, we show only the model fits, where we have used the distance of each model to convert the PMs from mas yr −1 to km s −1 so that isotropy is at 1; as we do not have the same data coverage (and, hence, bins) for the LOS and PM samples we cannot calculate these anisotropies directly, but the models given us an insight into what we cannot measure easily from the existing data. All panels are not constant, indicating that there is some anisotropy in the system and it varies through the cluster.
The right column shows, from top to bottom, the surface brightness profile, the projected M/L, and a proxy for the projected spherical anisotropy profile. The surface brightness profile is a key part of the model and so we are able to compare the model fit to the data, and it is a very good fit overall.
We cannot measure the mass directly -indeed, this is one of the main motivators for carrying out these models -so the M/L profile has no data points against which to compare. The M/L is clearly not constant through the cluster. It is ∼ 1.5 M /L at the centre, falls slightly in the intermediate regions and then increases to ∼ 10 M /L in the outer parts. This is consistent with a cluster that has some mass segregation, whereby high-mass stars tend to be more centrally concentrated than low-mass stars. The difference in mass between high-mass stars and low-mass stars is only around a factor of a few, 10 at most, but the low-mass stars are far more numerous than the high-mass stars, so it is the low-mass stars that dominate the mass budget. However, the high-mass stars are many hundreds or thousands times brighter than the low-mass stars, so it is the high-mass stars that dominate the light budget, despite their lower numbers. This interplay typically gives an M/L profile much like what we see here. The dashed line shows the global M/L of ∼3.0 for the cluster calculated as the total mass divided by the total light, which is consistent with stellar population synthesis estimates for clusters of its metallicity. This over-predicts the true M/L in the inner regions and under-predicts the true M/L in the outer regions.
Our best estimate for the total mass of the cluster is 2.38 +0.37 −0.30 × 10 6 M . This value is higher than the values reported by other authors. McLaughlin & van der Marel (2005) determined a mass in the range of (1.45 +0.28 −0.25 − 1.86 +0.33 −0.28 ) × 10 6 M based on M/L values of stellar population fits. Baumgardt & Hilker (2018) report a mass of 1.2 × 10 6 M based on a comparison between isotropic N-body models of globular clusters with observed LOS velocity measurements. We advise caution in interpreting these global measurements. The profiles are best constrained where there is a lot of kinematic  Figure 9. Best-fitting model profiles. From left to right then top to bottom: radial PM dispersion, tangential PM dispersion, LOS velocity dispersion, tangential PM / radial PM anisotropy, radial PM / LOS velocity anisotropy, tangential PM / LOS velocity anisotropy; surface brightness, projected M/L, proxy for projected spherical anisotropy. In all panels, the solid black line is the median of the fits, the darker shaded regions span the 15.9 to 84.1 percentiles (approximately the 1-sigma confidence region) and the lighter shaded regions span the 2.5 to 97.5 percentiles (approximately the 2-sigma confidence region). The coloured data points show data from this work and from literature sources to which the fits were performed as shown. In the anisotropy panels, the dotted lines highlight isotropy. In the mass-to-light panel, the dashed line shows the global M/L for the cluster calculated as the total mass over the total light. Overall, the fits to the data points are good. The cluster is best fit by a variable M/L, consistent with a cluster with some mass segregation, and a variable anisotropy profile, consistent with a cluster that is more relaxed in the central regions than the outer parts due to their different relaxation times. data (the intermediate radii) and less well constrained where the data is sparse (the innermost and outermost regions). In particular, the outer regions of the mass profile are driven by the outermost datapoint in the (Baumgardt & Hilker 2018) LOS profile, where there is no corresponding PM data at all. This is reflected in the large scatter we see in the outer regions of the fits.
The middle panels of Figure 9 show how the dispersions measured in two orthogonal directions compare in a 3-dimensional system. The lower panel in the right column effectively shows the 3-dimensional (spherical) anisotropy. Recall, in the JAM models, each component of the stellar density MGE is assigned an anisotropy . We calculate a luminosity-weighted anisotropy profile by multiplying the stellar density weights by the anisotropy and calculating the corresponding MGE profile. We divide this by the stellar density (surface brightness) MGE to get a proxy for the anisotropy profile. Here isotropy is at 0. The slight negative bias at the centre indicates some mild (spherical) tangential anisotropy there, although the models are also consistent with isotropy. The positive bias in the outer parts indicates (spherical) radial anisotropy there. This trend is consistent with theoretical expectations for cluster evolution and with cluster simulations. Clusters are expected to develop with some radial anisotropy, but they become more isotropic as they relax (Baumgardt & Makino 2003;Vesperini et al. 2014;Tiongco et al. 2016). As relaxation times are shorter at the centres, the centres relax and reach isotropy first. For NGC 6441, Harris (2010) reports a half-mass relaxation time of log( hm /yr) = 9.09, while the core relaxation time is significantly shorter with log( core /yr) = 7.93.
In addition, we made use of the mass profile determined with our model to recalculate the relaxation times. To do so, we used the equations from Djorgovski (1993) with the modification of Harris (2010). We took the core radius value of core = 0.13 arcsec from Harris (2010) (converted to core = 0.48 pc using our dynamical distance estimate), and assumed an average stellar mass of 1 3 M . The other parameters are based on our model. With our total mass estimate of 2.38 +0.37 −0.30 × 10 6 M , we find a half-mass radius of 8.0 +1.0 −0.7 pc. We determine a core density of (1.35 +0.19 −0.13 ) ×10 5 M pc −3 . The resulting values for the relaxation times are log( hm /yr) = 10.16 +0.11 −0.08 in the core, and log( core /yr) = 7.84 +0.02 −0.02 at the half-mass radius. While our values for the relaxation time in the core are similar to the values in the literature, our value for the half-mass relaxation time is larger as the previously reported values mainly because of the larger estimate for the half-mass radius. Due to the variable M/L ratio and our high total mass estimate, this value is higher then the half-light radius, which was used as a proxy for the half-mass radius in previous studies.
Although we advise caution in interpreting these global values (see above), we can see that the relaxation time at the half-mass radius is longer than in the core. This difference in relaxation times could explain why we observe isotropy in the centre, but anisotropy in the outer regions of the cluster. Similar astrometric findings for other clusters are discussed in Watkins et al. (2015a). Figure 10 shows the distribution of distance (top panel) and IMBH mass (bottom panel) estimates at the end of the MCMC run. To construct these plots, we selected 10 000 points in total, 100 walkers from each of 100 steps at 50 step intervals. The distance is nicely constrained, and approximately Gaussian overall. To determine a distance estimate and uncertainty, we calculate the median and 15.9 and 84.1 percentiles of the distribution to obtain distance from the values of 11.6 kpc from Harris (2010) and of 11.83±0.14 kpc from Baumgardt et al. (2019). The putative IMBH mass is not at all well constrained as evidenced by the flat distribution in the lower panel of Figure 10. With this data we cannot conclusively state whether or not an IMBH is present as models with a fairly massive BH and models with sub-solar-mass (effectively no) BH are equally likely, although models with very massive BHs are ruled out. The best we can do here is a place an upper limit on the mass of a possible IMBH at BH < 1.32×10 4 M .

Context of our IMBH mass limit
There are several scaling relations for black-hole masses in galactic bulges that are valid over a wide mass range. We extrapolate them to the mass of NGC 6441 and compare the predicted black hole mass with our upper limit.
Using our estimate for the mass of the cluster of 2.38 +0.37 −0.30 × 10 6 M , the extrapolation of the relation between the Bulge and BH masses of Schutte et al. (2019) gives an expected black hole mass of 1.17 +0.23 −0.18 × 10 3 M , significantly smaller than our upper limit. If we use the relation between velocity dispersion and BH mass of Gebhardt et al. (2000), our central velocity dispersion of (19.1 ± 2.0) km s −1 translates in an IMBH mass of 1.8 +0.37 −0.30 ×10 4 M which is quite similar to our upper limit.
Although these estimates are consistent with our upper limit, they are extrapolated from measurements of black-hole masses in a different (higher) mass regime than that of IMBHs in GCs. A more similar mass range is found in nuclear star clusters, in which the mass ratio between stellar cluster mass and black hole is typically BH / NSC ≈ 0.25 with a scatter of 2 (Nguyen et al. 2018;Greene et al. 2019). This is much higher than the upper limit of BH / NGC 6441 ≤ 0.0054 we obtain for NGC 6441. Theoretical predictions for IMBH masses in globular clusters have quite a range range between 0.001 and 0.01 of the cluster mass (Miller & Hamilton 2002;Portegies Zwart & McMillan 2002;Giersz et al. 2015), which is in agreement with our measured upper limit.
Finally Tremou et al. (2018) report an upper limit of M BH < 2270 M , based on the (missing) radio signature of an accreting IMBH in NGC 6441.
We can conclude that a larger sample of stars with precise PMs in the central region is necessary to observationally reach the lower end of the predicted black hole mass range and to put a final answer on the existence/non-existence of an IMBH in this cluster.

CONCLUSIONS
We determined PMs of around 1400 stars in the inner 15 arcseconds of the globular cluster NGC 6441 by combining space-based and ground-based position measurements taken 15 years apart. Because of the high astrometric precision of both epochs and the long timebase between them, our PM measurements reach a precision of 30 µas yr −1 for bright stars ( S < 13.5 or F555W < 17.5). Similar precisions are reached in pure HST-based studies (see Bellini et al. 2014) while the Gaia DR2 PMs have a higher uncertainty of around 0.1 mas yr −1 for stars with = 17 mag in less crowded fields (Lindegren et al. 2018). This proves the potential of combined ground-based AO and space-based astrometry with a long temporal baseline.
With the PM data we were able to determine the velocity dispersion profile of evolved stars in core of the cluster. In the innermost arcsecond we measure a velocity dispersion of (0.316 ± 0.034) mas yr −1 which corresponds to (19.1 ± 2.0) km s −1 assuming a distance of 12.74 kpc.
Using our PM measurements, we searched for signatures of a potential IMBH. Although we find one fast-moving object, whose projected trajectory is compatible with being ejected from the core because of the interaction with an IMBH, a statistical analysis of the PMs in our field does not show any signs of the presence of an IMBH in the centre of NGC 6441. A complete 5D picture (including radial velocity measurements) and deeper observations of the core are still needed to clearly confirm or rule out any hypothesis.
We used Jeans models to fit a combination of our newly-obtained kinematic data and existing kinematic catalogues. From the best-fit models we could determine the underlying physical properties of the cluster: the global M/L of the cluster is ∼ 3.0 M /L , but it varies from ∼ 1.5 M /L in the core of the cluster to ∼ 10.0 M /L in the outskirts. This is consistent with mass segregation where bright, high-mass stars are more centrally concentrated than low-mass stars. In the core we do not observe significant anisotropy, however the outer parts of the cluster show some radial anisotropy. By combining LOS velocity measurements with our PM measurements we obtain a dynamical distance of = 12.74 +0.16 −0.15 kpc. The models include a possible IMBH in the centre of the cluster. Our results are compatible with both the existence and non-existence of such a black hole, and we can only place an upper limit of BH < 1.32 × 10 4 M on the mass of the black hole. This value is about one order of magnitude larger than the mass predicted by extrapolating the relation between BH and bulge masses, but consistent with other predictions for GCs.
In future studies, deeper observations of the cluster would be beneficial as the number of stars in our kinematic sample is clearly limited by the number of detectable stars in the NACO exposures.
In the second half of this decade, instruments at extremely large telescopes, such as ELT MICADO, will allow PM measurements with an even higher precision and be able to finally answer the question of whether IMBHs are present in the cores of globular clusters.

APPENDIX A: GEOMETRIC DISTORTION CORRECTION
A set of geometric distortion corrections (GDCs) for NACO has already been published by Plewa et al. (2015Plewa et al. ( , 2018. However, the authors themselves noted that the distortion of the NACO detector is not stable over time, but shows abrupt changes most likely linked to instrument interventions. Furthermore, we do not know the exact definition of the PSF centre used to determine the literature GDC and there is a degeneracy between the PSF definition and the GDC. For these reasons, we decided to independently solve for the geometric distortion (GD) of the NACO detector using our data in order to achieve the best astrometric precision possible.
To determine the GDC, we tried both an autocalibration approach and the use of the ACS/HRC catalogue as a distortion-free reference. For the scientific analysis, we relied on the calibration based on the external ACS/HRC catalogue, as the resulting geometric distortion model led to a much better correction with fewer residual distortions.

A1.1 The reference catalogue
The HST ACS/HRC observations of NGC 6441 have a very high precision and a very reliable distortion correction that reaches an accuracy < 0.01 HRC pixel (see the Instrument Science Report Anderson & King 2004). In comparison with the uncorrected NACO catalogues, they can be considered effectively distortion free. However, we have to take into account the 15-year long time baseline between the HST and the NACO observations. The velocity dispersion in the cluster centre of around 18 km s −1 at a distance of 12.74 kpc (our dynamical distance estimate) leads to an rms displacement of around 0.4 NACO pixel. While this effect is purely statistical and is averaged out when measurements of multiple stars are combined, it still leads to a decreased precision of the GDC and can mask smaller GD effects.
To overcome the limitations caused by the stellar motions, we made use of a PM catalogue for the core of NGC 6441 created using the HRC exposures from 2003 and WFC3/UVIS observations from 2014 and that will be the subject of a future paper (Bellini et al., in preparation). Suffice here to say that the data reduction and proper-motion computation of this catalogue closely followed the prescriptions given in great details in Bellini et al. (2014) and Bellini et al. (2018b). The number of well-measured stars in the core of the cluster is relatively small due to the larger pixel scale of the WFC3/UVIS channel (40 mas pixel −1 ). On the other hand, we only use bright, well-measured stars to determine the GDC anyway.
The PMs in this catalogue were used to propagate the HRC 2003 positions to the epoch of the NACO observations. We only included stars in the reference catalogue with total PM error PM < 0.07 mas yr −1 , which corresponds to an error in the displacement of 0.08 NACO pixel. Furthermore, we restricted the selection of reference stars using different quality criteria. In the end, around 1600 stars were available as reference stars.

A1.2 Matching individual frames on reference catalogue
As a first step, we matched the astrometric catalogues containing uncorrected NACO positions onto the reference catalogue using linear transformations. The parameters of the linear transformations were determined using a least-squares fit. The cut-off radius for matching stars was initially 3.5 NACO pixels and it was progressively decreased down to 1 NACO pixel during the iterative process.
To avoid absorbing the linear distortion terms (the so-called skew) already at this level, we used 4-parameter transformations, which contain a shift in x and y direction, a rotation and a change of scale.
The residuals of the linear-transformation fit now contained both the individual measurement errors, and also the geometric distortion of the NACO images. We divided the detector in 10 × 10 quadratic bins and collected the residuals from all matched NACO catalogues. We then calculated the 3 -clipped median of the residuals in each bin containing at least 20 residuals. The maps of the binned residuals can be seen in Figure A1 and provide a first indicator of how the geometric distortion correction will look like.

A1.3 Fit of a 2D polynomial model
We modelled the binned residuals using the following 2D third order standard polynomials: = 1 + 2 + 3 2 + 4 + 5 2 + 6 3 + 7 2 + 8 2 + 9 3 = 1 + 2 + 3 2 + 4 + 5 2 + 6 3 + 7 2 + 8 2 + 9 3 (A1) We used a least-square fit to determine the coefficients. The coefficients 1 and 2 were set to zero to lower the degrees of freedom of the fit and to enforce that, at the centre of the detector, our GDC will lead to the same -scale as the detector and the corrected and raw -axis will be aligned. However, as we cannot assume that and axis have the same scale/are perpendicular, we can not set the parameters 1 and 2 to zero. A higher order of the polynomials did not lead to a better fit of the data. Also, the use of a different polynomial base (e.g. Zernike polynomials) had no significant influence on the result.

A1.4 Iterative Process
To avoid that the linear transformations between the single NACO catalogues and the reference frame are biased by uncorrected geometric distortion, we repeated the procedures described in an iterative process. After the polynomial coefficients have been determined, we applied 75% of the corrections to the raw NACO catalogues and repeated the full process (described in A1.2 and A1.3) with the new corrected coordinates until the polynomial coefficients converge (in our case after 200 iterations).

A2 Determination of the GDC with an autocalibration approach
Initially, we planned to solve for the GDC using an autocalibration approach in which the distortion-free reference frame is obtained by combining multiple different pointings. As the position measurements of each star are based on measurements at different parts of the detector, the effect of the GD randomises and therefore is averaged out. This is a well-proven technique and has been used to calibrate the GD of the HST (Anderson & King 2003;Bellini & Bedin 2009), but also for various ground-based studies (Anderson et al. 2006;Bellini & Bedin 2010;Libralato et al. 2014Libralato et al. , 2015. We refer to these publications for a detailed description of the iterative process. In comparison to the GDC obtained with an external reference (see section above) our autocalibration result showed significant differences (see Figure A2): we were unable to determine the linear distortion (the so-called skew) terms, which caused a rms difference of 0.68 NACO pixel, but also the nonlinear differences had an rms of 0.28 NACO pixel.
By comparing the NACO master frames with independent HST catalogues of NGC 6441 (based on WFC3/UVIS and ACS/WFC observations), we could verify that the differences between the distortion corrections indeed are caused by uncorrected GD in the autocalibration reference frame.
It can be easily understood why we were not able to determine the linear distortion terms using the autocalibration: the linear skew terms are the same over the whole field of the detector. Even if stars are measured on different detector positions, their position measurements are affected by the same skew and therefore the reference frame also has the same skew. This degeneracy is lifted if an instrument with multiple detectors is used, as is the case for the papers cited above, or if there are pointings with different orientation on sky.
The nonlinear differences are most likely caused by a relatively low number of stars in comparison with e.g., the autocalibration of the HST or ground-based wide-field instruments.
The failure of the autocalibration approach with our 2018 NACO dataset highlights the importance of the choice of the dither + rotation pattern. Especially if no external reference with sufficient astrometric precision is available (as will possibly be the case for future instrument at ELTs), a particular attention on this is required during the preparation of future observations. Table B1 contains the values of the PM-based velocity dispersions computed in our work.  0.298 ± 0.016 0.273 ± 0.020 0.329 ± 0.024 11.75 0.287 ± 0.016 0.285 ± 0.023 0.289 ± 0.023 Table B1. Results for the proper motion dispersion profiles in the inner region of NGC 6441.