Detecting Globular Cluster Tidal Extensions with Bayesian Inference: I. Analysis of $\omega$ Centauri with Gaia EDR3

We base our analysis on high-quality proper motions and photometry provided by the Gaia mission, and in particular EDR3 (Gaia Collaboration et al., 2016; Gaia Collaboration et al., 2020). We began by retrieving all stars within a five degree radius of $\omega$ Cen from the Gaia archive¹¹1https://gea.esac.esa.int/archive/ using Table Access Protocol (TAP) facilities. With our adopted distance to the cluster of 5.2 kpc (Harris, 1996; Soltis et al., 2021), this corresponds to a physical radius of $\sim$450 pc. We first applied the several adjustments to the data as suggested by the EDR3 documentation. Specifically, these are the zero-point correction in the parallax using the code provided in Lindegren et al. (2020) and the correction of the $G$-band magnitude and the corrected flux excess factor as provided in Riello et al. (2020). We then applied a series of cuts to the resulting catalogue, namely:

• •

  stars with $(G_{BP}-G_{RP})>1.6$ mag are excised as this part of the colour magnitude diagram (CMD) primarily contains field dwarfs;

• •

  due to overestimation of the $G_{BP}$ flux causing sources to appear to blue, we removed all stars $G_{BP}>20.3$ mag (see Section 8.1 in Riello et al., 2020);

• •

  stars whose corrected BP and RP excess flux is greater than 3 times the associated uncertainty (see Section 9.4 in Riello et al., 2020) are removed as this typically relates to poor photometry;

• •

  stars with poor astrometric data are removed, using the ‘Re-normalised Unit Weight Error’ , $u_{r}>1.4$. This value has been shown to provide the optimal selection of stars with ‘good’ astrometric solutions (Lindegren et al., 2018, 2020)

• •

  stars are removed that have a well-measured parallax that places them within 3 kpc of the Sun: $(\varpi-3\sigma_{\varpi})>0.3$ mas. This removes all nearby stars that contaminate the $\omega$ Cen sightline.

These selections pruned the original sample of 2.5 million stars to 1.7 million stars, a removal of 30 per cent. We then transformed all positions and proper motions onto a tangential coordinate projection about the centre of $\omega$ Cen using the equations (2) in Gaia Collaboration et al. (2018c) - ($\xi$,$\eta$) now denote positions (for $\alpha$ and $\delta$, respectively) and $\mu^{*}_{\xi}$ and $\mu_{\eta}$ denote the associated proper motions. In addition, we also corrected the proper motion of all stars in the sample for the solar reflex motion at the adopted distance of 5.2 kpc using the python gala package (Price-Whelan, 2017), which uses the solar motion $(12.9,245.6,7.78)$ km s^-1 (Drimmel & Poggio, 2018). Any stars that do not have a proper motion measurement in EDR3 are also removed at this stage.

For the subsequent analysis, we estimate EDR3 photometric uncertainties from the flux zero-point equations given in Riello et al. (2020). Furthermore, the photometry is de-reddened using the Schlafly & Finkbeiner (2011) dust maps and the colour-dependent extinction equations given in Gaia Collaboration et al. (2018b). Given the large area of the sky spanned by our dataset and the highly variable extinction across this region (see Fig. 1), the de-reddening correction has been performed on a star-by-star basis, and the resultant de-reddened photometry is denoted as $G_{0}$, $G_{BP,0}$ and $G_{RP,0}$.

After the initial cuts on the data were performed, we proceeded to remove further contaminants from our sample through additional photometric selection. Specifically, we wanted to retain only those stars that were consistent with the well-defined sequence for $\omega$ Cen in colour-magnitude space (Fig 2, left). Because incompleteness due to crowding affects the quality and depth of the photometry in central regions of the cluster, we defined a ‘fiducial sample’ of stars that lie within a cluster-centric radius of 10 to 50 arcmin. The outer radius of this sample is just beyond the nominal King tidal radius of 70.25 pc presented in de Boer et al. (2019), which corresponds to 46.4 arcmin at our adopted distance. The proper motion of $\omega$ Cen has been derived to be ($-3.25,-6.75$) mas yr^-1 by Vasiliev & Baumgardt (2021) using Gaia EDR3. We refine the fiducial sample by removing stars with proper motions that differ by more than 1 mas yr^-1 of this value when transformed to the solar reflex motion-corrected frame: $(\mu^{*}_{\alpha},\mu_{\delta})$=($3.08,-3.58$) mas yr^-1. To guide our photometric selection, we fit a PARSEC²²2http://stev.oapd.inaf.it/cmd isochrone (Bressan et al., 2012; Marigo et al., 2017) with [Fe/H]$=-1.53$ dex, [$\alpha$/Fe]$=+0.2$ and an age of 12 Gyr in the Gaia EDR3 bandpasses to the fidicual sample. Although it is established that multiple stellar populations exist within $\omega$ Cen which results in significant colour width on the RGB, these are believed to be largely concentrated within the central 10 arcmin (e.g., Bellini et al., 2009) and we assume that they have negligible effect in the outermost regions that we focus on here. We assign each star in the fiducial sample a pseudo-colour value, $w$, which is defined as the absolute difference between the $(G_{BP}-G_{RP})_{0}$ colour of the star and that of the isochrone at the corresponding $G_{0}$, which we denote as $(G_{BP}-G_{RP})_{ISO}$, divided by the star’s uncertainty in colour $\sigma_{(G_{BP}-G_{RP})0}$ (see also Kuzma et al. 2018). That is:

We determine the distribution of $w$ as function of $G_{0}$ to identify the best $w$ range that describes the cluster sequence, assuming a half-normal distribution. At a given $G_{0}$, we record the standard deviation of the distribution ($\sigma_{w}$) and adopt $2\sigma_{w}$ to capture the bulk of the stellar sequence. The function $\sigma_{2}(G_{0})$ then represents how $2\sigma_{w}$ varies as a function of $G_{0}$. Due to systematic effects, $\sigma_{2}(G_{0})$ increases significantly at the bright end (i.e. as $G_{0}$ and the photometric uncertainties decrease) so we employ an exponential fit to $\sigma_{2}(G_{0})$ to allow for this effect. This keeps the red giant branch (RGB) sufficiently covered while selecting almost all stars on the main sequence and the turn-off regions, minimizing the contamination from MW field stars. Finally, with $\sigma_{2}(G_{0})$ defined, we calculate $w$ for stars within the entire (5^∘ radius) sample, and keep only stars that satisfy the condition $w\leq\sigma_{2}(G_{0})$. This removes a further 1.2 million stars from the sample, and yields a final sample of $5\times 10^{5}$ stars that we feed to our model. The fitted isochrone. the $w$ range and the boundary of the $\sigma_{2}(G_{0})$ cut are shown in Fig. 2 (right). We note that this selection excludes blue horizontal branch (BHB) stars and RR-Lyrae (RRL) stars at this point in the analysis but we will revisit those populations in a later section.

Our Bayesian technique has measured the posterior distributions for 15 different parameters, and the solutions are presented in Table 1 (the posterior distributions are presented as on-line material). The non-zero normalisation constants derived for $f_{cl+ex}$ and $f_{cl}$ indicate that our analysis has recovered the cluster as well as detected an extended extra-tidal feature surrounding it.

In Fig. 3, we present the proper motion distribution of our sample, highlighting those stars with an $\omega$ Cen membership probability of $\geq 0.5$. The cluster clearly stands out relative to the dominant MW field distribution. We fit the solar reflex motion corrected proper motion of $\omega$ Cen to be (${\mu^{*}_{\xi}},\mu_{\eta})=(3.055\pm 0.002,-3.624\pm 0.002)$ mas yr^-1. This corresponds to (${\mu^{*}_{\xi}},\mu_{\eta})=(-3.21,-6.74)$ in the observed frame and is in excellent agreement with the recently determined EDR3 values from Vasiliev & Baumgardt (2021). We find that $f_{cl+ex}\approx 0.10$, which indicates that ${\sim}10$ per cent of our photometrically-selected stars belong to the cluster + extended structure components on the basis of their proper motions and spatial positions, while the remaining ${\sim}90$ per cent belong to the MW field. This high level of MW contamination is not surprising given the low Galactic latitude of $\omega$ Cen and the large area on the sky that we have analysed. Nonetheless, the membership probability distribution presented in Fig. 4 demonstrates that our method is effective in distinguishing between cluster and contaminant stars at both small and large cluster-centric radii, with a couple thousand high probability members lying beyond the King tidal radius.

We unambiguously detect an extra-tidal component to $\omega$ Cen in the form of compelling tidal tails which span the $\sim 10^{\circ}$ on the sky that we have analysed. The normalization factor between the cluster and the extended structure is $f_{cl}=0.910\pm 0.002$, indicating that ${\sim}91$ per cent of the stars we assign to the cluster reside within the King tidal radius, while ${\sim}9$ per cent lie in the tidal extension. However, we stress that this is not representative of the overall fractional mass or luminosity in the extended structure component since photometric completeness issues at small and large radii, as well as the stellar mass function of stars in the different components, are complicating issues (see Section 4.2). The position angle of the extended component is well defined, $\theta_{ex}=122.88^{\circ}\pm 2.14^{\circ}$ (East of North), indicating a preferred orientation. Notably, this is not in the same direction as the direction of the MW field gradient, which is $\theta_{MW}=147.85^{\circ}\pm 0.35^{\circ}$ (East of North).

The morphology of the extended component is visualised in Fig. 5 which shows the surface density distribution of stars with $P_{mem}\geq 0.5$. Stars are assigned to 3 arcmin $\crossproduct$ 3 arcmin bins on the sky and subsequently smoothed with a Gaussian of width 12 arcmin. We calculated the mean bin density of stars lying outside 1.2 times the Jacobi radius 161.7 pc or $106$ arcmin (Balbinot & Gieles, 2018) and the associated spread, hence removing the influence of the cluster itself in our statistics. Each bin is then displayed as the number of standard deviations above the aforementioned mean bin density. For comparison, we also show the distribution of all 500,000 stars, created under the same conditions, that passed our initial selection criteria prior to the Bayesian analysis. The dominant gradient across the field seen here has been effectively removed by our technique to reveal the underlying structure. 5 also shows the model fits to the extended structure (quadrupole) and MW contamination (linear), demonstrating they describe the data well.

The tidal tails extend across the full extent of the analysed area, corresponding to a physical diameter of $\sim 900$ pc. We have explored greater radial extents surrounding the cluster but find no significant detection of the debris beyond our initial search radius of 5 degrees; however, on such large scales, our assumption of a single value for the cluster proper motion will begin to break down. The tidal tails we have uncovered exhibit the same overall morphology as those seen in other recent studies of $\omega$ Cen (Ibata et al., 2019a; Sollima, 2020), as well as those hinted at in earlier work (Marconi et al., 2014), but have been recovered here through a completely independent method and one which allows membership probabilities to be associated to individual stars. Fig. 5 also shows the orbital path of $\omega$ Cen calculated with Galpy³³3http://github.com/jobovy/galpy, using the MWpotential2014 from Bovy (2015) and the $\omega$ Cen positional and velocity values from (Baumgardt et al., 2019). There is a reasonably good correspondence between path of the orbit and the direction of the extended structure, which is to be expected as tidal extensions are typically found to closely follow the orbit (e.g., Montuori et al., 2007).

Fig. 6 shows the radial fall-off of the extended component derived from counting stars along the direction of the tails (on-axis) and perpendicular to it (off-axis). For constructing these profiles, we have used only $P_{mem}\geq 0.5$ stars that are within the $2\sigma$ density detection contour of Fig. 5. The on-axis profiles were calculated using stars that lie within a 45 $\deg$ wedge, centered on $\omega$ Cen, of the position angle defined by $\theta_{ex}$. The off-axis profiles included those stars within a 45 $\deg$ wedge at an angle perpendicular to $\theta_{ex}$. Within these wedges, we counted the number of stars in concentric annuli out to 5 degrees, with the size of the annuli increasing at large radii to combat small number statistics. Due to issues with incompleteness, we exclude stars in inner 20 arcmin of $\omega$ Cen (see also de Boer et al., 2019) and instead represent the behaviour of the radial profile within this radius using the azimuthally-averaged V-band surface brightness profile of Trager et al. (1995). The surface brightness profile and star count profiles are stitched together by scaling the datapoints in the 20-40 arcmin range, allowing us to construct a radial profile that spans $\sim$ 20 magnitudes in surface brightness. Poisson uncertainties are shown for the star counts.

The on and off axis profiles exhibit very similar behaviour to a radial distance of $\sim$ 40 arcmin ($\sim 60$ pc), indicating a predominantly spherical morphology out to the King tidal radius. Beyond this radius, the stellar distribution becomes increasingly elongated along the axis of the tails, as manifest by both a higher surface density of stars at a given circular radius as well as detections further out. Also shown in Fig. 6 are power-law fits to both the on and off-axis profiles for $P_{mem}\geq 0.5$ stars lying beyond the King tidal radius. For the on-axis profile, we find that the density follows a $\gamma_{\rm on}=-3.40\pm 0.20$ decline, and a sharper decline for the off-axis profile, $\gamma_{\rm off}=-5.22\pm 0.26$. These values are consistent with the power-laws seen in clusters with tidal tails originating from mass-loss (e.g., Carballo-Bello et al., 2014; Küpper et al., 2015).

The fraction of the present-day mass in the tidal extensions can be calculated by integrating the radial profile in Fig. 6 using the LIMEPY code (Gieles & Zocchi, 2015). For this calculation, we assume that mass follows light and that the mass function of stars is a constant in these peripheral parts. To facilitate comparison with other work, we use the Wilson truncation radius, found to be 70.6 $\pm 1.2$ arcmin (consistent with de Boer et al. 2019), as our fiducial radius. For reference, this is $\sim 1.5$ times the King tidal radius. Taking the average of the on and off-axis spherical integrations, we find that the amount of mass in the radial range from the Wilson truncation radius to the edge of our field, is a mere $0.1$ per cent of the total mass in the system. Interestingly, this is considerably less than the $\sim 1$ per cent mass fractions contained in the stellar envelopes seen around NGC 1851 and NGC 7089, which have been calculated in an identical manner (Kuzma et al., 2018). It is also far less than the mass fractions of 3–50 per cent seen in the tidal tails of some lower mass MW globular clusters, such as Palomar 5 (Odenkirchen et al., 2003), NGC 5466 (Grillmair & Johnson, 2006) and M92 (Thomas et al., 2020). However, it is worth noting that our estimate is based on only the stream stars which lie within a five degree radius of the cluster, and hence is likely to be a lower limit on the fractional mass in the extensions.

We also explore how the ellipticity and position angle of $\omega$ Cen vary as a function of radius from the main body into the tails. To do this, we fit elliptical isophotes to the 2D surface density distribution displayed in Fig. 5 using the fitting routine in the python package, photutils. This package performs the fitting based on the iterative method introduced by Jedrzejewski (1987) and in our analysis we choose to hold the centre fixed. The resulting profiles are presented in Fig. 7 and show that $\omega$ Cen remains roughly spherical to a radius of about 1.5 degrees, which corresponds to just inside the Jacobi radius, before becoming increasingly elliptical at larger distances. The position angle converges on $\theta_{ex}$ at approximately 1.5 degrees as well.

Prior to our work, the ellipticity of $\omega$ Cen had only been explored out to 30 arcmin, which is well within the King tidal radius. Within this radius, Calamida et al. (2020) have shown that the cluster becomes increasingly elliptical as the radius decreases, peaking at a value of $\sim 0.16$ at 8 arcmin. Our measurements at radii $\lesssim 20$ arcmin are unreliable due to incompleteness issues with Gaia but in the range of $20-30$ arcmin we are in excellent agreement with Calamida et al. (2020). Future data releases of Gaia will be better suited to deal with the highly crowded regions of GCs and will allow homogeneous study of the ellipticity and position angle of $\omega$ Cen from its inner regions to furthest extent.

In our analysis, we have only considered cluster stars that lie on the upper main sequence, main sequence turn-off and RGB regions of the CMD. However, Fig. 2 demonstrates that $\omega$ Cen possesses a prominent horizontal branch as well, composed of BHB and RRL stars. To examine whether the distribution of these populations is also consistent with the tidal extensions, we revisit our final sample of $5\times 10^{5}$ stars described in Section 2.2 and select those objects with $14\lesssim G_{0}\lesssim 16$ and $-0.3\lesssim(G_{bp}-G_{rp})_{0}\lesssim 0.4$. We also require that these stars have proper motions that lie within 3 sigma of the cluster value listed in Table 1. We verified that those objects remaining have parallaxes consistent with that of $\omega$ Cen. Fig. 5 shows the locations of these candidate BHB stars superposed on the surface density distribution of main sequence and RGB stars. The majority of the BHB candidates are confined to the main body of the cluster, which is to be expected given that it is the lower mass (i.e., main sequence) stars that are preferentially populating the tidal extensions (Balbinot & Gieles, 2018). However, those BHB candidates that do lie outside the cluster region follow the tidal extensions rather closely.

To explore the RRL star distribution, we relied on the tables vari_rrlyrae and vari_classifier_result from the Gaia DR2 analysis of Holl et al. (2018). We selected sources that lie within a five degree radius of $\omega$ Cen and, as for the BHBs, that have proper motions that lie within 3 sigma of the cluster value. As shown in Fig. 5, the RRL candidates are also highly clustered in the central regions but the three objects which lie beyond the Jacobi radius show a good alignment with the tidal features. For these three objects, we calculate an time-averaged dereddened magnitude of $G_{0}=14.2\pm 0.02$. Assuming the absolute magnitude of RRL (type AB) to be $M_{G}=0.64\pm 0.25$ (Iorio & Belokurov, 2019), we calculate the distance to these stars to be $5.2\pm 0.6$ kpc, which is entirely consistent with the measured distance of $\omega$ Cen (Harris, 1996; Soltis et al., 2021). The small number of candidate RRL stars that we have found coincident with the $\omega$ Cen tidal debris is in agreement with the earlier work of Fernández-Trincado et al. (2015), who failed to find any strong RRL candidates in a search area of 50 sq. degrees around the cluster. Of the small number of potential candidates they identified beyond the tidal radius, we have used EDR3 to confirm that none of these stars have proper motions consistent with $\omega$ Cen membership.

Our analysis is based on a mixture model approach for spatial positions and proper motions, applied to a sample of stars that have been pre-selected on photometric properties. In this section, we examine what radial velocity (RV) data exist for our sample and whether these measurements are consistent with our membership assignments.

The RV of $\omega$ Cen is well-determined and large (234 km s^-1 Baumgardt et al., 2019), meaning that there is excellent contrast with the surrounding field population. Gaia DR2 (and EDR3) includes mean RVs for 7.2 million stars brighter than $G\approx 13$ and with $T_{\rm eff}$ in the range $\sim 3550-6900$ K (Gaia Collaboration et al., 2016; Gaia Collaboration et al., 2020). We explored which of our sample of $5\times 10^{5}$ stars in the vicinity of $\omega$ Cen have Gaia RVs. Of the 48 stars that do, we find that 13 of these stars have $P_{mem}\approx 1$ while the remaining 35 have $P_{mem}\approx 0$. Reassuringly, the $P_{mem}\approx 1$ stars all have RVs which lie within 10 km s^-1 of the mean $\omega$ Cen velocity whereas the $P_{mem}\approx 0$ all have RVs $\lesssim 170$ km s^-1, further confirming they are non-members. All the stars with Gaia RVs lie within 30 arcmin of the cluster centre.

We also searched other public spectroscopic survey datasets for RVs for our sample and found that 295 of our stars have RVs in APOGEE DR16 (Jönsson et al., 2020). Of these, we find 291 have $P_{mem}\approx 1$ and RVs consistent with cluster membership but that only two of these lie beyond 30 arcmin radius. These numbers are unsurprising given that the APOGEE-2 observations come from a targetted program on $\omega$ Cen but it is still encouraging that so many are returned as high probability members by our method.

There have been very few dedicated spectroscopic searches for stars in the far outer regions of $\omega$ Cen stars thus far. Da Costa & Coleman (2008) and Da Costa (2012) (hereafter DC12) used the AAT 2dF to search for candidate stars on the lower red giant branch over the radial range 20-60 arcmin, i.e. from roughly half to just beyond the King tidal radius. We cross-matched the probable members and probable non-members from DC12 with our catalogue. We find that our sample contains 102 out of the 160 DC12 probable members, of which we confirm 99 as high probability members ($P_{mem}>0.5$). Those not present in our sample, as well as those not recovered as high probability members, are all located within the inner 30 arcmin of the cluster, where crowding affects the astrometry and photometry and where the multiple population signature causes some stars to be missed by our photometric selection. Interestingly, our technique also identifies 11 ( 30 per cent) of DC12’s probable non-members as stars with $P_{mem}>0.5$. These turn out to be stars that have RVs consistent with that of $\omega$ Cen but that DC12 suspected were contaminants on the basis of their line-strengths. Our analysis indicates that these 11 stars are indeed members of the cluster. All of the cross-matched member stars lie within the tidal radius.

Finally, Sollima et al. (2009) conducted a search using VLT/FLAMES for $\omega$ Cen stars at large radii. Cross-matching our catalogue with their RV members, we find 90 stars in common, all of which have $P_{mem}\approx 1$ and lie within 35 arcmin of the cluster centre.

The fact that our algorithm can independently recover known RV members with (very) high probability is reassuring but unfortunately the comparisons we have been able to make are quite limited due to the fact that there are so few known RV members beyond the tidal radius. Indeed, the majority of targets we have identified as members of the tidal features are fainter than the typical depth of existing spectroscopic surveys and the heavy contamination along this sightline means that previous targetted programmes have been inefficient in identifying RV members at large radii. Our sample of $\approx 4000$ high probability candidates outside 30 arcmin provides an excellent sample for future spectroscopic follow-up, which will enable confirmation of membership as well as measurements of RVs and chemistry.