Discovery of Astrometric Accelerations by Dark Companions in the Globular Cluster ω Centauri

We present results from the search for astrometric accelerations of stars in ω Centauri using 13 yr of regularly scheduled Hubble Space Telescope WFC3/UVIS calibration observations in the cluster core. The high-precision astrometry of ∼160,000 sources was searched for significant deviations from linear proper motion. This led to the discovery of four cluster members and one foreground field star with compelling acceleration patterns. We interpreted them as the result of the gravitational pull by an invisible companion and determined preliminary Keplerian orbit parameters, including the companion’s mass. For the cluster members, our analysis suggests periods ranging from 8.8 to 19+ yr and dark companions in the mass range of ∼0.7 to ∼1.4M ☉. At least one companion could exceed the upper mass boundary of white dwarfs and can be classified as a neutron star candidate.


INTRODUCTION
Since the beginning of measuring proper motions it is implied that individual stars move along straight lines.This is still a good approximation for most sources even now in the Gaia era.If the motion indicates curvilinear pattern, that would change the proper motion by astrometric acceleration, usually measured in mas yr −2 .Stronger accelerations may allow to invoke the Keplerian formalism (e.g.Jia et al. 2019).Thus, the detection of astrometric accelerations near the radio source Sagittarius A * was instrumental to eliminate all other alternatives to the supermassive black hole (Ghez et al. 2000;GRAVITY Collaboration et al. 2018).This discovery and additional astrometric accelerating sources Corresponding author: Imants Platais, Johannes Sahlmann imants@jhu.edu,Johannes.Sahlmann@ext.esa.int in the same area are examples of ground-based astrometry of binary motion where one component is dark.
In 1990-1993 the Hipparcos space astrometry mission regularly observed multiple times ∼118,000 pre-selected stars.For the first time in history regular astrometric observations were obtained outside the Earth's atmosphere and covering the entire celestial sphere.This effort resulted in 2622 acceleration candidates (Lindegren et al. 1997) but none with a genuine invisible companion.
The third data release (DR3) of the Gaia space astrometry mission used data acquired between July 2014 and May 2017 and reported 338,215 acceleration solutions and 165,500 orbital solutions (Halbwachs et al. 2023).A few illustrative cases of compact and invisible companions are presented by the Gaia consortium in Gaia Collaboration et al. (2023a).To 'weigh' an invisible companion one must determine the mass of the primary source and it is desirable to obtain epoch radial velocities.This recipe made it possible to discover Gaia BH1 and Gaia BH2, two dormant black holes (BH) in intermediate-period binaries (El-Badry et al. 2023a,b).An opportunity to enhance the accuracy of proper motions and to probe very long-period accelerations is the combination of the Hipparcos and Gaia catalogs (e.g.Brandt 2021).
Although the Hubble Space Telescope (HST ) is not an astrometric observatory, it has provided spectacular science results such as the detection of the very small tangential velocity of the Andromeda galaxy at 17.0 km s −1 (van der Marel et al. 2012), establishing the upper limit of a putative, intermediate-mass black hole in ω Cen (van der Marel & Anderson 2010), precision parallaxes (±45 µas) for long-period Cepheids in the Milky Way with the spatial scanning technique (Riess et al. 2018), and the detection of astrometric microlensing by a black hole (Sahu et al. 2022).
These breakthroughs gave us confidence that measurements of astrometric accelerations are feasible with the prime HST WFC3/UVIS camera.One configuration that can generate a measurable acceleration in a few years is a double star with a faint or dark secondary component.In order to maximize the number of surveyed sources and therefore the odds of detection, we focused on dense stellar system such as the core of globular clusters.Thanks to the foresight of the HST WFC3 team, a multi-cycle calibration program targeting the center of ω Cen (NGC 5139) was set in motion in 2010.So far, a total of ∼200 well time-spread frames have been accumulated, spanning over a total of 13 years.In certain aspects these observations may surpass the available Gaia catalogues due to the 2-dimensional pointed observations, the longer time span, and the slightly better spatial resolution which is critical in very crowded stellar environments such as the cores of globular clusters.However, the Gaia limitations due to source blending and onboard resources discussed in (Gaia Collaboration et al. 2016, Sect. 6.6) were partially mitigated by a special observing mode (e.g.Sahlmann et al. 2016) and a focused product release in the ω Cen region was published on 10 October 2023 (Gaia Collaboration et al. 2023b), see Sect.6.3.
A similar HST program (GO-12911, PI L.R. Bedin) targeted the nearest (∼1.7 kpc) Galactic globular cluster M 4 (NGC 6121 Bedin et al. 2013).However, the timespan of these observations is only one year with an expectation that this would be sufficient to detect the astrometric 'wobble'.As of this writing, no accelerating source in M 4 has been reported.The combined analysis of proper motions from HST and Gaia DR3 claim a central dark mass of 800 M ⊙ in M 4 (Vitral et al. 2023).

HST /WFC3 OBSERVATIONS OF ω CEN
The observations of ω Cen were taken with the fourthgeneration imaging camera WFC3/UVIS on-board HST installed during Servicing Mission 4 in May 2009, with the angular field of view of 2. ′ 73 × 2. ′ 73.The first observation in calibration program CAL-12094 (PI J. Kim Quijano) was taken on 2009-07-15 to monitor the flatfield uniformity through the filter F606W at the adopted center of ω Cen (RA=201.• 69283 Dec=−47.• 47905).The following calibration programs CAL-11911 (PI E. Sabbi), CAL-12094 (PI L. Petro), CAL-12339 (PI E. Sabbi) and the successive calibration programs 12353,12714,13100,13570,14031,14393,15000,15593,15733,16413,16588 (all by PI V. Kozhurina-Platais) are addressing a variety of specific aspects of the WFC3 instrumental calibration, such as photometric stability, geometric distortions, and their stability (for the latter aspect, see also Sahlmann et al. 2019).All of these programs used F606W besides some other filters and relatively short exposure times of 35-40 or 60 sec and a variety of roll angles from −170 • to 180 • .During the first five years nearly 50% of frames had a large 40 ′′ dither.From the view-point of precision astrometry such dithers could be detrimental due to the partial overlap.The large variety of roll-angles produced a lower number of total detections outside the radius of r=80 ′′ from the adopted center.At larger radius and certain roll angles some targets cannot be observed.This is visible in Figure 1.Altogether, for this study we selected a total of 175 images spanning 13 years with the last observation on 2022-01-11.
The distribution of epochs over each calendar year is not even.There are two intervals without observations: April 29 to June 3 (35 days) and September 6 to December 12 (97 days).These 'holes' are due to the unschedulability of ω Cen at those times.Some of the other missing observations are due to failure of guidestar acquisition.In some cases, such periodic pattern may spawn spurious one-year signals in the periodogram analysis (Section 6).

DATA REDUCTION
The chosen flc-type images from the HST Mikulski Archive for Space Telescopes (MAST) are a standard output of the HST calibration pipeline CALWFC3 with bias and dark subtraction, and flat-fielded on original pixels.The flc.fits images are corrected for the charge-transfer efficiency (CTE) in the HST pipeline (Anderson et al. 2021).The effect of CTE tend to increase with time and can affect the point-spread function (PSF).This can have an impact on the measurements of stellar centroids and flux.Therefore, the pixelbased corrections is a crucial step towards deriving high- hst1pass produces the catalog of X&Y positions, the measured flux provided in the instrumental magnitude and qfit, the quality parameter of the PSF fit.We limited our sample to measurements with 0.0<qfit≤0.1.The typical distribution of qfit as a function of instrumental magnitude can be found in (Figure 4; Platais et al. 2015), which shows two easily-separable populations of artifacts (with larger qfit) and genuine measurements.Our goal is quality, not completeness, and the applied constraint only selects the very-best measurements.The derived high-precision positions are corrected for geometric distortions (Bellini et al. 2011) and hst1pass also simultaneously provided celestial coordinate RA & Dec.The WFC3/UVIS geometric distortions are known to be stable over time (Kozhurina-Platais & Anderson 2015).Nevertheless, the still-ongoing distortion monitoring is one of the goals in the multi-cycle WFC3 calibration program.

Astrometric registration
The output of hst1pass provides a geometricallycorrect gnomonic projection for a tiny part of the celestial sphere.We used Gaia Early Data Release 3 (EDR3) (Gaia Collaboration et al. 2021) to provide the reference frame.Preliminary tests showed that EDR3 source quality is atypical in the center part of ω Cen.If we limit the EDR3 sources to the ones with positional error less than 1 mas, then nearly 50% of them are missing proper motion and there is a cut-off at G ∼18 mag.Clearly, this is a consequence of the crowding and limitations of Gaia onboard resources (Gaia Collaboration et al. 2016;Cantat-Gaudin et al. 2023).
Similar to Bellini et al. (2017) we selected the seedframe ibc301qrq flc taken on 2010-01-14.If we translate the Gaia celestial coordinates to the epoch 2010.12 and then convert into gnomonic coordinates and align them to the measured pixel coordinates of the seed frame, the rms of the resulting residuals is ∼4.4 mas.There are 440 common stars for this transformation.In contrast, any pair of WFC3/UVIS frames with close epochs generates an rms of ∼0.7 mas.Such disparity may imply that the EDR3 positions near the center of ω Cen have substantially lower accuracy then expected.For this reason we decided to use the EDR3 positions only for aligning the seed-frame to the RA and Dec axes and to derive the WFC3/UVIS pixel scale in the F606W filter.
To align the seed frame it had to be rotated by 149.
• 751 around the location of star 227418 (Table 7 Bellini et al. 2017).At this stage we decided to use the Bellini et al. (2017) numbering system to designate our sources since it covers the same area.Also, the availability of celestial coordinates in both data-sets make it easy to crossidentify the sources.Our estimate of the WFC3/UVIS pixel scale is 39.770 mas and is not corrected for the velocity aberration.As noted by Bellini et al. (2011);Anderson (2007) the effect of velocity aberration changes only the pixel scale.In this sense, our pixel scale is not absolute.
Next, all frames were converted into pixel-based coordinates of the seed-frame using a least-squares minimization.The number of common stars in the magnitude range 15.4 < F606W < 17.5 span from 1093 to 2434 depending on the dither pattern and the HST roll angle.If the pixel coordinates are free of geometric distortions, then a linear 3-term polynomial should be sufficient to align any pair of frames.Our analysis of the residuals yielded the following results: • The rms of the residuals is correlated with the epoch of observations.Early epochs produce low rms at ∼0.015 pix whereas the last epochs have values as high as ∼0.22 pix.The significant inter-nal velocity dispersion near the center of ω Cen at 0.75 mas yr −1 (Anderson & van der Marel 2010) is the main reason this increase.Such patterns are not only detrimental to the astrometric accuracy but may also impact the prospects of finding accelerations.
• Initial solutions using only three linear terms revealed a small (up to 0.05 pix) but persistent quadratic pattern in the residuals, mainly along the Y -axis.Unfortunately, we do not know whether the seed-frame itself is free from quadratic markings.There are a handful of cases with nearly perfect distribution of residuals but not necessarily at the same HST roll angle.
• About a dozen frames show atypical box-shaped distributions of residuals.Since such frames appear close in time to typical frames, we suspect that sometime the guiding of HST was not perfect, especially for some short exposures used in this calibration program.
The solution to these issues is a three-pronged approach: • Subtract from all X&Y pixel-positions the effect of proper motion relative to the epoch of the seedframe.We subtracted the effect of the known relative proper motion from each X&Y position at every epoch relative to the epoch of the seed-frame.
In practice, the known proper motion is multiplied by the time difference in years and the result is subtracted from X&Y .As the result, all positions have zero proper motion up to its formal precision.
Then, these positions are averaged and the mean is subtracted.In essence, our approach levels the positional accuracy over all epochs.This can be seen in Figure 2 in the even scatter of residuals in time.
At each epoch we calculate the amount of this effect and convert it into pixels using exact scale at 40 mas per pixel and then rotate around the adopted zeropoint at X c =2073.19 and Y c =2292.38 pix using the HST position angle P A V 3. Since the proper motion is strictly vectorial, it cannot affect the inferred acceleration.The opposite, however, is true: acceleration can bias the inferred proper motion.
• The modified pixel coordinates X m , Y m are then translated into the system of the seed-frame positions by a least-squares regression using a polynomial model with 6 terms (linear and quadratic) for both axes (e.g., Bellini et al. 2014).The output are epoch positions X f , Y f free of proper-motion effects to the level of their accuracy.For sources of particular interest, e.g. with likely acceleration signals, we also subtracted a small local offset in both coordinates by using all sources within ∼ 50 pixels around the target and at all epochs.Then the averages of these offsets are applied to all epochs.
• A total of 17 epochs are excluded from the final analysis, the recalculation of proper motion, and the search for acceleration.It is most likely that at these epochs HST experienced jitter in its pointing stability due to the degrading gyro (Anderson & Sabbi 2018).There is no sharp divide between the normal jitter (∼ 11 mas) and the increased ones.We note that such HST calibration frames may not be science grade but still can be used for engineering purpose.Hence they are in the MAST.

LONG-PERIOD ACCELERATIONS
The new astrometric catalog includes 162,604 sources down to m F 606W ≃ 21.5.We examined the parameter space so that it is optimal for the detection of accelerations.The suitable limits appear to be at m F 606W 16.1 − 20.0 mag and no less then 120 epochs evenly spread over 13 years.There are 22,466 such sources.Each set of X f ,Y f then is used to calculate small additional proper motion and acceleration which can be approximated by the quadratic term in each axis, X 2 and Y 2 .Thus, the output of the least-squares minimization yields two formal accelerations along each axis, named ξ −2 X and ξ −2 Y (see Table 1).The crucial parameter is the ratio the acceleration term and its formal error, designated as Q X and Q Y .This ratio (ignoring the sign) should be at least ∼5 in one or both axis (see Q X and Q Y in Table 1 and Figure 2).With this criterion we find four clear cases of astrometric acceleration discovered in a globular cluster.
Using a smaller threshold of 3-5 on the Q-parameter generates ∼400 candidates, but these are heavily polluted by impostors.A typical example is an offset or unusual spread at some narrow bunch of epochs in the same HST program.Generally, our search is limited by the known noise-floor of geometric-distortion corrections at ∼0.3 mas (Bellini et al. 2011).
In addition, the long-period acceleration of source 233697 was detected in the independent periodogram analysis of the astrometric time-series, see Section 6.Among the astrometric accelerations this source has the shortest period at 8.8 yrs.Due to the complete orbit and more, the quadratic term will diminish or may even flip its sign.We inspected the F606W photometric timeseries and found them to be stable with robust scatter estimates (Lindegren et al. 2012) of between 14 and 42 mmag, without noticeable correlations with the acceleration signals.

MASSES OF VISIBLE COMPONENTS FROM STELLAR EVOLUTIONARY MODELS
Three of the cluster members with high acceleration are clearly main sequence stars located just below the turn-off of ω Cen, see Fig. 4, and the fourth one (ID=233697) is on the sub-giant branch.
Table 1.Astrometric binaries discovered in this work: astrometry at the Gaia DR1 reference epoch J2015.0 (from Bellini et al. 2017) and photometry of the visible components and acceleration parameters.mV and mI are HST WFC3 Vega magnitudes in F606W and F814W filters; proper motions µX , µY , and their errors σµ X , σµ Y in mas yr −1 ; N frame number of frames (epochs) followed by accelerations and their significance.The last column indicates our designation of these newly-discovered binaries.

ID
12.9 14.7 HSToC-4 AB -0.0098 -0.0274 2.8 7.7 HSToC-5 AB In the absence of any information other than the photometry, their mass can only be estimated via theoretical isochrones (Serenelli et al. 2021).For instance, one can simply select the isochrone with the most appropriate value of age and metallicity, then shift it to the appropriate distance modulus and foreground reddening and then read the stellar mass from the isochrone that is closest to the observed star in the color-magnitude diagram (CMD).
In the case of ω Cen this procedure is complicated due to the wide range of metallicity (see e.g.Frinchaboy et al. 2002).Indeed, if we take a 12-Gyr isochrone with the mean metallicity of ω Cen, namely [Fe/H]= −1.53, and relocate it to the distance of 5.24 kpc (Soltis et al. 2021) and redden it by E(B − V ) = 0.132 mag (cf.Bono et al. 2019), then the high-acceleration stars clearly fall to the red, meaning that they are probably more metal rich.
In order to create a more realistic model for these stars, we simulated a cluster composed of single stars only, with an age of 12 Gyr and using the entire distri- bution of [Fe/H] derived by Frinchaboy et al. (2002).We note that their metallicity distribution function is based on photometry with Washington M , T 2 and DDO51 filters, but that it was confirmed by the high-resolution spectroscopy from APOGEE (Mészáros et al. 2021).This simulation contains a tail of metal-rich stars extending far enough to the red and cover our four accelerating stars (see left panel in Fig. 5).To estimate the masses of the four visible components with high acceleration, we proceed as follows: For every star in the simulation with color c and magnitude m, a weight w i is assigned as where (c 0 , m 0 ) is the color-magnitude position of the visible components, and σ = 0.02 mag is the typical photometric error for both color and magnitude.These weights reflect the likelihood that the colormagnitude location of the observed stars derive from the theoretically-predicted location plus a Gaussian distribution of photometric errors.These weights are used together with the masses and metallicities of the simulated stars, to build the probability density functions (PDF) for both mass and metallicity.These functions are depicted in the central and right panels of Fig. 5. From these PDFs, we then derive the median and 68% confidence (1σ-equivalent) intervals for both mass and metallicity.Their values are listed in Table 2. Table 2 initially presents the values obtained by assuming the standard reference to the distance of ω Cen (5.24±0.11kpc; Soltis et al. 2021) based on Gaia EDR3 parallaxes.The error in distance modulus in this case is of just 0.004 mag, which is one fifth of the error we already assume in our derivation of the masses and metallicities, and hence not able to significantly affect our estimations.However, there are alternative values to ω Cen distance that are worth exploring.The review of GC distances by Baumgardt & Vasiliev (2021), for instance, provides a 5.426 ± 0.047 kpc distance; these authors argue that the Gaia parallaxes have small-scale correlated errors and systematic biases but this is not well yet understood in EDR3, particularly for ω Cen.The bottom part of Table 2 presents the mass and metallicity estimates that are obtained with this distance value.As it can be appreciated, changes in the mass values are smaller than 0.01M ⊙ .
En passant, we note that the location of our highacceleration stars on the simulated CMD indicate that the metallicities of these stars are higher than the mean value of ω Cen [Fe/H]= −1.53 (see Table 2).

PERIODOGRAMS OF ASTROMETRIC TIME-SERIES
To search for deviations from the linear astrometric motion expected for a single source, we processed the two-dimensional astrometric timeseries as follows.For every source we subtracted the average position from the timeseries and fitted a four-parameter linear model corresponding to two positional offsets and two proper motions.This analysis is therefore independent of the amount of proper-motion displacement subtracted in the pre-processing and analysis discussed in Sect.3.1.We omitted the parallax term because we do not expect to measure the parallax signal for cluster members, due to the inherently relative nature of our astrometric measurements.
We then computed the periodogram of the fit residuals using the kepmodel package (Delisle & Ségransan 2022).The period and false-alarm probability (FAP) of the highest periodogram peak were retained and are shown in Figure 6.The FAP indicates the probability that the peak is caused by random noise, i.e. a small FAP indicates the presence of periodicity in the timeseries.The lower and upper boundaries of the periodicity search were 10 days and 20 000 days, respectively.
Most sources exhibit insignificant periodogram peaks with large FAP, and we performed detailed analysis for all sources with small FAP.A few sources show significant signals with periods of one year or half a year.These are likely foreground sources for which measurable parallax offsets are present.The list of sources with FAP < 10 −10 and periods around half a year is 420482,414263,396157,421049,421381,427775,428247,428253,428571,434388,412665,421379,421520,427880,396159,421074.The list of sources with FAP < 10 −10 and periods around one year is 247164, 276969, 109223, 87685, 387712, 116672, 107882, 312186.As a by-product of our analysis, we have therefore identified 24 sources that probably are foreground sources and not members of the Omega Cen cluster.
The pile-up at long periods is caused by the upper boundary of the period search.Some of these peaks have small FAP but visual inspection of the time-series revealed systematic offsets and/or few measurements in the astrometric timeseries.The same applies to the few short-period (P < 100 day) sources and to the small-FAP source at P ∼ 10 yr visible in Figure 6.
Finally, we identified several sources with long-period periodogram peaks and small FAP whose origin is genuine Keplerian motion.

Constraints on the orbital parameters and masses of the companions
For the sources identified in the previous section on the basis of their small periodogram FAP, we attempted to derive constraints on the possible orbital parameters.One source (143023) turned out to be a foreground binary and is discussed in Appendix A. The periodograms of the astrometric timeseries residuals for the remaining four sources are shown in Fig. 7.All of them except 233697 show a dominant peak at roughly the observation timespan with a large tail towards longer periods.This is typical for astrometric data that cover only a small or moderate portion of the astrometric orbit.This is also evident when inspecting the timeseries in Fig. 2.
For source 233697, whose timeseries is shown in Fig. 3, the periodogram peak is at ∼ 3230 days (8.85 yrs) with an observation timespan of ∼4730 days (12.96 yrs), indicating that more than one orbital revolution has been observed.
We analyse the data of these sources under the assumption that the non-standard motion is caused by a single dark companion in a Keplerian orbit about the visible source.In three out of four cases, the observation time-span is too short to derive tight constraints on the orbital parameters and we therefore attempt to determine the range of possible configurations and their astrophysical interpretations.
For every source we use kepmodel to fit the data with the initial four-parameter model plus a 7-parameter Keplerian model with initial parameters estimated from the periodogram.The orbit-fitting results are therefore independent of the amount of proper-motion displacement subtracted in the pre-processing (Sect.3.1).We do not impose priors on any of the parameters but we set an upper limit of 0.3 to the eccentricity for 143023 and 317591.For the remaining sources that limit is 0.95.This is to avoid high-eccentricity solutions which often are preferred in orbit-fitting results using incomplete orbital coverage or otherwise sub-optimal sampling.We assigned a uniform uncertainty of 0.3 mas to every individual astrometric measurement, which is motivated by the distortion-correction noise floor.This choice is validated by the typical value of an additional jitter term, which we included in the model.The model fitting is always performed on the individual-frame data, hence no binning is applied.For visual display, however, we compute normal points as averages of individual frames grouped in time.
We fit the combined model with a standard maximumlikelihood algorithm and report a few of the best-fit parameters in Table 3, which has some caveats: • The algorithm tends to converge towards the smallest period compatible with the data set.However, the true period of the system may be much longer which influences the estimated parameters of companion.
• For two sources we applied the e < 0.3 constraint.
The true eccentricities may be larger than 0.3, which can imply larger semi-major axe and, therefore, companion mass estimates.Table 3 also shows the solution parameters for larger eccentricities.For 143023 a larger eccentricity results in a larger semimajor axis and longer period, but the companion mass is nearly unchanged.For 317591 a larger eccentricity results in a longer period and larger companion mass, but the semimajor axis is nearly unchanged.
• When inferring the companion mass, we assumed one single dark companion.This assumption may not be justified, which we discuss below.The alternative scenarios include a pair of quasi-dark companions or a companion object that is not dark but contributes significant light.In the latter case, the mass of that companion has to be even larger than what we inferred.
The uncertainties of orbital period P , eccentricity e, and semi-major axis of the primary a 1 (Table 3) correspond to the formal uncertainties returned by the minimisation routine and are indicative only.The companion mass M 2 was estimated from the orbital parameters and an average of the primary-mass estimates M 1 in Table 2 1 , with uncertainties propagated using Monte-Carlo resampling.
Since several of these orbits are not fully covered by our data and both proper motions and orbital parameters are free parameters in the fitted model, it can happen that proper motion is mistaken for orbital motion and vice versa.This is however not the case for our solutions.We verified that the residual proper motion determined as part of the orbital fit is comparable to or smaller than its uncertainties.This is important to safeguard against potential instabilities in the fitting routine where a spurious proper motion and long-period Keplerian offsets compensate each other.
Although the systematic and formal uncertainties are significant, we can categorize the four binaries into two groups: 233697 and 317591 are probably systems with similar component masses, whereas the solutions for 212028 and 290133 indicate companions that appear to be significantly more massive than the primary star, but yet might be much dimmer.
The secondary masses in Table 3 were derived under the assumption that the companion is dark, i.e. the positions of the system's photocentre and the primary coincide.The semimajor axis of the photocentre motion a photo that we determine directly is then equal to the semimajor axis of the primary's barycentric motion a 1 = a photo .If the companion contributes light, the photocentre motion gets diluted and a 1 > a photo which implies an even larger mass of the companion.The scaling between a binary's photocentre and barycentre orbit size is determined by the fractional mass f = M 2 /(M 1 +M 2 ) of the components and their fractional luminosity where āphoto is the photocentre orbit size in linear units and ārel is the semimajor axis of the relative orbit in linear units.Using Kepler's third law in SI units this can be written as and therefore (4) From Eq. 4 it is clear that to maintain the observed a photo constant between the dark companion limit (where L 2 = 0 and a 1 = a photo ) and a luminous companion with L 2 > 0, the companion mass M 2 has to increase.
With the data available to us we cannot exclude that the companions contribute light, but if they do it implies that they are fainter and at the same time more massive than the primary.
Two configurations that could satisfy these constraints for HSToC-2 and HSToC-4, which have the largest estimated companion masses, are (a) the companion is a tight pair of white dwarfs (see Section 8) or (b) the companion is a tight pair of main-sequence stars.We simulated the latter configuration by generating synthetic triple systems that have the same astrometric signature in terms of period but are composed of a primary with properties from Table 2 and a twin-binary companion of the same metallicity with individual masses between ∼0.4-0.8M ⊙ .Using the isochrones discussed in Section 5 and Eq. 2 we found that the corresponding photocenter orbit size is largest when the total companion mass is ∼1 M ⊙ , but that it remains 0.2 − 0.25 mas smaller than the measured orbit size.This scenario therefore seems excluded for HSToC-4 and unlikely for HSToC-2.

Orbital motion of 290133/HSToC-4
The source 290133 stands out as having the largest estimated companion mass.We show the orbital motion of 290133 in the plane of the sky in Figure 8 and as a function of time in Figure 9.For the sake of better visualisation only, we computed averages over time-bins of ∼100 days (normal points).We stress that these figures show only one possible realisation of the orbital parameters, the actual orbit may look different.However, these visualisations are useful to ascertain the Keplerian orbital motion.
In Figure 9 it can be appreciated that the excess signal after subtracting parallax and proper motion has a standard deviation of 0.51 mas (both axes, computed on normal points) which is reduced to 0.25 mas when including the Keplerian model.The equivalent figures for the other sources are shown in Figures 15, 17, and 16 of the Appendix.

Gaia catalogues in ω Cen
In the Gaia DR3 source catalog there are no counterparts to the binary sources listed in Table 1 when we impose a maximum separation of 5 ′′ , a Gaia magnitude m G < 20 and magnitude difference |m V − m G | < 2, and the presence of a Gaia parallax measurement.
On 10 October 2023, however, Gaia published a Focused Product Release (FPR) that includes an catalogue of additional sources in ω Cen (Gaia Collaboration et al. 2023b).In that catalogue (Gaia archive table gaiafpr.crowdedfield source, epoch J2017.5),we found unique positional matches to all five sources within a radius of 30 mas.We did not correct for proper motion between the DR1 and FPR epochs.The details of these matched sources are given in Table 4.They have small magnitude differences |m V − m G | < 0.3 mag and the next closest positional matches are at separations ≳ 0.3 ′′ .We briefly explore the properties of these Gaia data here.
Figure 10 shows the FPR 'astrometric excess noise' (AEN) parameter for our binaries in comparison with roughly comparable sources in that catalogue.This parameter corresponds to a jitter term that is added quadratically to the individual Gaia measurement uncertainties to account for the residual scatter in the Table 3.Preliminary estimates of some Keplerian parameters and the masses of the companions.These were obtained under the assumption of a single and dark companion and with a prior on eccentricity for sources HSToC-1/143023 and HSToC-5/317591.Formal uncertainties from the minimisation routine are given only for orientation.M1 is derived from Table 2.The foreground source HSToC-1/143023 is discussed in Appendix A. The second part of the table shows indicative results for HSToC-1/143023 and HSToC-5/317591 when relaxing the eccentricity prior and allowing for values of ∼0.5 and ∼0.8.The KRV is the estimated radial-velocity amplitude.single-source astrometric solution (Lindegren et al. 2021).Without access to the Gaia astrometric timeseries it is a powerful metric to identify and sometimes even characterise astrometric binaries (e.g.Kiefer et al. 2019;Belokurov et al. 2020;Gandhi et al. 2022).Whereas HSToC-1/143023, HSToC-3/233697, and HSToC-5/317591 exhibit values close to the mode of the distribution, HSToC-2/212028 and HSToC-4/290133 show elevated levels of Gaia excess noise.These latter sources are also the ones with the larger orbital signatures (a 1 ≃ 0.9 mas) inferred from the HST data compared to HSToC-3 and HSToC-5 with a 1 ≃ 0.5 mas.This is a strong indication that the Gaia FPR is sensitive to orbital motion with amplitudes of ∼1 mas in ω Cen.The foreground source HSToC-1 does, however, not match that pattern since we estimated an orbit size of ∼0.9 mas, but its AEN is unremarkable.This may be due to the long (and uncertain) period of this binary compared to the FPR timespan of 5 years, which reduces the probability that orbital acceleration shows up as elevated AEN.
It is evident in Figure 10 that there are many sources with AEN values similar or even much larger than HSToC-3 and HSToC-5.A fraction of those are likely of instrumental origin, e.g.source blending, but another fraction could indicate acceleration signals.Using the Gaia FPR astrometric excess noise may therefore be a promising avenue for the discovery of new astrometric binaries in ω Cen, especially in regions not covered by the HST observations.

MUSE OBSERVATIONS
The sources 212028 and 290133 are covered by the footprint of the mosaic observed with MUSE (Bacon et al. 2010) as part of the survey of Galactic globular clusters presented in Kamann et al. (2018).The quality cuts routinely performed prior to searches for binary stars (see Giesers et al. 2019) result in sets of six radial velocity (RV) measurements for 212028 and five RV measurements for 290133, with typical uncertainties of > 10 km s −1 for the former and 5 − 10 km s −1 for the latter.None of the stars shows signs of variability in the MUSE data.However, given the relative large uncertainties compared to the expected velocity semiamplitudes (Table 3), this does not seem unexpected.
We note that the spatial resolution of the MUSE data, with a typical FWHM of ∼ 0.8 arcsec, is lower compared to the HST data by a factor of ∼8.Therefore, even though the stars are bright enough to perform accurate RV measurements with MUSE, crowding limits our abilities to do so.Both sources have similarly bright stars located at about 0.5× the FWHM and significantly brighter neighbours at distances of ∼ 1 ′′ .Future observations with the MUSE narrow field mode will be able to overcome this problem.

DISCUSSION AND FUTURE STUDIES
The largest Milky Way globular cluster ω Cen now is the first cluster with detected astrometric accelerations.Whether or not the relatively small range of the deduced masses of the invisible components is typical for Milky Way globular clusters remains to be answered by future studies.However, we can get some cues from another rare populations such as millisecond pulsars (MSP) and X-ray sources.Curiously, the first bona fide MSPs in ω Cen were discovered only in 2020 (Dai et al. 2020) followed by observations with the MeerKAT radio telescope and discovery of 13 more MSPs (Chen et al. 2023).Four of them are in our field (Figure 11).
Since MSPs and neutron stars share the same physical phenomenon observed in different wavelengths, we inspected the known binary systems anywhere in the sky with well-measured component masses (Table 1 Özel & Freire 2016).These binary systems yield fairly-tight mass limits of 1.38 ± 0.05M ⊙.Two of our targets, 212028 and 290133, fit very well to these mass limits and thus provide a strong argument that their invisible component could be a neutron star.In fact, the estimated properties of these two systems are very similar in luminosity, acceleration, component masses, and Keplerian parameters such as the semi-major axis.
An alternative option for the nature of these companions would be a tight pair of white dwarfs, but such a scenario may not be effective in the cores of globular clusters.The presence of MSPs in our field itself is auspicious to the neutron star scenario.
Old stellar systems such as globular clusters contain low-luminosity X-ray sources.The most comprehensive rework of archival ω Cen observations with the Chandra X-ray Observatory Advanced CCD Imaging Spectrometer provide the best sample of these sources (Henleywillis et al. 2018).There are 14 X-ray sources within the boundaries of our field.None of them is close to the detected astrometric binaries.In globular clusters the typical X-ray sources are cataclysmic variables (Belloni & Rivera 2021).They are relatively tight binaries composed of a white dwarf or a neutron star and a main sequence star as donor.Our field has four cataclysmic variables but their potential semi-major axis would be totally imperceptible.Future studies can expand on our work in several ways.Including archival HST ACS/WFC images of ω Cen taken prior to WFC3/UVIS images and accumulating more data from the ongoing WFC3/UVIS calibration program will expand the time baseline of the observations and therefore lead to better constraints on the binaries' orbital and mass parameters.
Independently, one could experiment with relaxing the conservative cutoff value for the image quality parameter qfit to increase the number of available measurements per source, hence potentially increase sensitivity.
Access to the epoch positions X f ,Y f , free of the effect of proper motions using Bellini et al. (2017) data, and new calibrated magnitudes will be granted upon request to the authors and made available through the SciServer platform.
Further exploration of the Gaia FPR dataset of ω Cen has the potential to lead to the discovery of additional binaries in this cluster.

SUMMARY AND CONCLUSIONS
This study is trailblazing astrometric acceleration searches in a relatively distant globular cluster populated by roughly solar-and subsolar-mass sources.We demonstrated the viability of this project by obtaining high-precision astrometric timeseries for tens of thousands of sources, thereby improving the proper motions of more than 22,000 cluster members with an average precision of 0.011 mas yr −1 .We discovered five new binaries, one in the foreground and four cluster-members, and characterised their Keplerian motion.
Despite the large distance of ω Cen of over 5 kpc we were able to detect the first astrometric binaries.Those have very-small orbital semi-major axes below one 1 mas (Table 3).These binaries have periods longer than our survey timespan, except for HSToC-3, where we observed more than one Keplerian orbit allowing us to derive better constraints on the mass of the companion.Inspection of the Gaia FPR dataset of ω Cen independently supports the discovery of the larger-amplitude binaries HSToC-2 and HSToC-4, which have elevated levels of Gaia astrometric excess noise.
Our preliminary best-fit mass estimates for the invisible companions in ω Cen range from ∼0.7M ⊙ to ∼1.4M ⊙ and from ∼0.6M ⊙ to ∼2.1M ⊙ if we account for the associated 1-σ uncertainties.The true range of allowed companion masses is likely even larger, but there are essentially three likely scenarios of their nature: a single white dwarf 2 , a pair of white dwarfs, or a neutron star.Because of the uncertainties in our estimates and the Chandrasekhar limit on the maximum mass of white dwarfs at 1.4 M ⊙ , we are unable to pinpoint the most-likely scenario.
Our initial goal was to probe the existence of stellarmass black holes in binaries of ω Cen.Although we have not unambiguously identified a BH candidate, we have demonstrated that surveys like ours are sufficiently sensitive to discover them if they exist.A system akin to Gaia BH2 at the distance of ω Cen generates an astrometric signal with a semi-amplitude of ∼0.9 mas and a period of 3.5 years, which lies well within the sensitivity limit of our survey, cf.Table 3. Shorter period systems akin to Gaia BH1 would generate a signal of ∼0.25 mas and a period of 0.5 years, thus are more challenging to detect.I.P. thanks Andrea Bellini for sharing the astrophotometric catalog, Timothy Brandt for testing a target, Terrence Girard for the parallax-related code and general advising.J.S. thanks Jari Kajava for checking for X-ray and radio counterparts to our binaries.
This research makes use of the SciServer platform.The SciServer is a collaborative research environment for large-scale data-driven science.It is being developed at, and administered by, the Institute for Data Intensive Engineering and Science at Johns Hopkins University.The authors gratefully acknowledge grant support for HST program AR-16629, provided by NASA through grants from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract 5-26555.The HST data underlying this research are available at MAST: 10.17909/dwtd-kr81.
This research made use of pystrometry, an open source Python package for astrometry timeseries analysis (Sahlmann 2019).
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.The astrometry periodogram of source HSToC-1/143023 after fitting the four-parameter linear model reveals a strong signature at a one-year period which persists in the residuals of tentative orbit-model fits.This is a strong indication that this is a foreground source whose un-modelled parallax motion introduces the periodicity.We therefore computed the parallax factors using the formalism provided by van de Kamp (1981) and then fitted the standard five-parameter model.As shown in Fig. 12 this cancels the power at one year.We therefore modelled the astrometry of this source with models that include the parallax term.
To estimate the potential orbital parameters, we modelled the data with the combination of the standard fiveparameter model and a single Keplerian model with initial parameters compatible with the highest periodogram-peak.Since the orbital motion is only partially covered (Fig. 2) we can only determine a lower limit to the period and give indicative bounds to the possible mass range of the companion.Since orbital eccentricity is essentially unconstrained at this point, we limited its allow range to <0.3 (see Table 3).
One possible orbital solution with a period of ∼ 7000 days is shown in Fig. 13.However, the orbital period could be as short as the ∼ 4000 day timespan of the data or much longer.The orbital motion and fit residuals for this solution as a function of time are shown in Figure 14.
The determined parallax of the system is almost independent of the orbital parameters because the timescales are very different for such long-period systems.We determined a relative parallax of ϖ rel = 0.55 ± 0.12 mas.Since that is measured relative to cluster members, we need to account for the cluster distance of 5.24 kpc (Soltis et al. 2021) and obtain an absolute parallax of HSToC-1 of ϖ abs = 0.73 ± 0.12 mas.Adopting this latter parallax value, and assuming that the HSToC-1 primary is a metal-poor halo star with the same foreground extinction as ω Cen, its placement in theoretical isochrones provides a mass estimate in the range between 0.65 and 0.71 M ⊙ .Assuming a primary mass of 0.68 , the companion mass for this particular solution would be ∼0.15M⊙ , where we neglect any light contribution by the companion, which could dilute the photocentre motion and lead to an underestimation of the companion mass (Table 3).Belokurov, V., Penoyre, Z., Oh, S., et al. 2020, MNRAS, 496, 1922

Figure 1 .
Figure 1.Distribution of the locations of two stars on the WFC3/UVIS CCD detector over 175 epochs.Red points correspond to star 290133 with 174 detections, blue points correspond to star 212028 with 149 detections.Red points are more concentrated and that resulted in higher overall accuracy.The background are stars with mF 606W ≤ 18.5.

Figure 2 .
Figure 2. Four cases of long-period acceleration.Magenta lines show quadratic least-square fits.Green lines indicate zero offset.If a star lacks acceleration, all measured XY (black points) would cluster around a straight line in both axes.As an example of borderline acceleration is star 317591 in X-axis.This star would be rejected if not the prominent acceleration in the Y -axis.

Figure 3 .
Figure 3. Star 233697 has a low acceleration signal, indistinguishable from normal stars.Its Keplerian motion was detected in the periodogram (Figure 7).Green line indicates zero offset, although a careful inspection shows small outbalances within 1 mas.

Figure 4 .
Figure 4. Observational CMD of ω Cen in the HST photometric Vega system through WFC3/UVIS F606W and F814W filters (Bellini et al. 2017).A randomly-chosen subset of ∼10% of all stars is shown.Red circles show the locations of the astrometric binaries.The reddest of these is the field binary HSToC-1.

Figure 5 .
Figure 5. Left panel: a simulated CMD of ω Cen (black tiny dots), built to estimate the masses and metallicities of the four visible components with high acceleration (coloured dots).The central and right panels show the probability density functions for both mass and metallicity for these stars, as derived from the stars with similar color-magnitude position in the simulation (see text for details).

Figure 6 .
Figure 6.Density histogram of FAP as a function of peak period for the ∼160 000 sources for which the periodogram could be computed.The two vertical lines are at periods of 0.5 and 1 year.The black stars indicate the location of the five sources in Table1.

Figure 7 .
Figure 7. Residual periodogram of the four-parameter model for 212028 (top), 233697 (second), 290133 (third), and 317591 (bottom).The vertical grey dashed line shows the location of the highest peak, whose period and FAP are indicated in the panel title.The black dotted line corresponds to the observation timespan for this source.

Figure 8 .
Figure 8. Orbital motion of HSToC-4/290133 in the sky after subtracting the standard 5-parameter model.The grey curve corresponds to the best-fit orbital parameters.Uncertainties represented as error-bars were estimated from the dispersion of the data accumulated in the normal point.Dashed lines connect observed and computed locations.

Figure 9 .Figure 10 .
Figure 9. Orbital motion of HSToC-4/290133 after subtracting the standard model (top panels) as a function of time.The final residuals of the full-model fit are shown in the middle and bottom panels.The top and middle panel show normal points computed by combining individual frames, whereas the bottom panels show individual frame data.Left panels are in RA direction and right panels are in Declination.

Figure 11 .
Figure 11.Distribution of rare objects in ω Cen.The center of the gnomonic projection is the best-known kinematic center of ω Cen.Black dots show randomly-chosen 7% of all program stars.Red points are oue newly-discovered binaries.The one in the upper-right quadrant is the foreground field binary.Green points indicate millisecond pulsars and blue ones X-ray sources.

Figure 12 .
Figure 12.Residual periodogram of the four-(top) and five-parameter (bottom) models for HSToC-1/143023.The power at long period is hardly affected but including the parallax in the model makes the power around 365 days disappear.

Figure 13 .
Figure 13.Parallax and proper motion (top) and orbital motion (bottom) of HSToC-1 in the sky.

Figure 14 .
Figure 14.Orbital motion of HSToC-1 after subtracting the standard model (top panels).The standard deviation of 0.59 mas (both axes, computed on normal points) in the top panels is is reduced to 0.32 mas when including the Keplerian model.

Figure 15 .
Figure 15.Orbital motion of HSToC-2 after subtracting the standard model (top panels).The standard deviation of 0.38 mas (both axes, computed on normal points) in the top panels is is reduced to 0.27 mas when including the Keplerian model.

Figure 16 .
Figure 16.Orbital motion of HSToC-3 after subtracting the standard model (top panels).The standard deviation of 0.30 mas (both axes, computed on normal points) in the top panels is is reduced to 0.16 mas when including the Keplerian model.

Figure 17 .
Figure 17.Orbital motion of HSToC-5 after subtracting the standard model (top panels).The standard deviation of 0.36 mas (both axes, computed on normal points) in the top panels is is reduced to 0.31 mas when including the Keplerian model.Belloni, D., & Rivera, L. 2021, in The Golden Age of Cataclysmic Variables and Related Objects V, Vol.2-7, 13, doi: 10.22323/1.368.0013

Table 2 .
Astrophysical parameters of the visible components.Masses were estimated on the basis of the sources' positions in the CMD.

Table 4 .
Gaia FPR identifiers and relevant properties.The ρ column indicates the match distance and the AEN column lists the FPR 'astrometric excess noise'.