© The Authors 2025

1 Introduction

The origin and structure of Omega Centauri (ω Cen) is a highly debated topic. Being the most massive and luminous globular cluster (GC) in the Milky Way, it has been argued that ω Cen may not, strictly speaking, be a GC but rather the remnant core of a tidally disrupted, nucleated dwarf galaxy (Bekki & Freeman 2003; Wirth et al. 2020; Johnson et al. 2020). Evidence in favor of this hypothesis includes the observation of a retrograde orbit (Dinescu et al. 1999), which may be associated with the disruption of the so-called “Sequoia” (Myeong et al. 2019) or Gaia-Enceladus galaxies (Massari et al. 2019), the detection of a long tidal stream (Ibata et al. 2019), and a new formation mechanism for nucleated dwarfs that was recently proposed by Gray et al. (2024). In their model, the smallest nucleated dwarfs first have their star formation quenched by reionization before it is reignited again in a major-merger-driven starburst. The dense nucleus forms in this starburst, yielding a galaxy with two distinct stellar populations that have a large age gap. This observational feature has been seen in several massive Milky Way GCs, including ω Cen (Gray et al. 2024).

Recent numerical simulations from Wirth et al. (2020) suggest that if ω Cen is indeed an accreted nucleated dwarf then it can retain significant amounts of dark matter even as it tidally disrupts. By performing a Jeans analysis of stellar kinematics, Brown et al. (2019) and Evans et al. (2022) found that the presence of an extended dark component in ω Cen is statistically favored, and that the self-annihilation of dark matter could account for the observation of γ-ray signals. However, millisecond pulsars (MSPs) are known γ-ray emitters and can also account for the observed γ-ray emission (Reynoso-Cordova et al. 2021; Dai et al. 2023, 2020).

Alongside the possible presence of dark matter, ω Cen is interesting as a possible location for the formation of black hole seeds via runaway stellar collisions (e.g., Portegies Zwart et al. 2004) or the formation and death of a supermassive star (e.g., Chon & Omukai 2020). Several authors (Noyola et al. 2008; Noyola et al. 2010; Jalali et al. 2012; Baumgardt 2017) found evidence of a ~3−5 × 10⁴ M_⊙ intermediate-mass black hole (IMBH) at the center of ω Cen from the consideration of stellar kinematic data (using either Jeans models or N-body simulations). However, Van der Marel & Anderson (2010), via a Jeans-based analysis, derived a 3σ upper bound on the IMBH mass of 1.8 × 10⁴ M_⊙, in tension with the previously cited masses. Baumgardt et al. (2019) also found no evidence of an IMBH by comparing N-body simulations to data for ω Cen (stellar velocity dispersion, surface brightness, and proper motion data), favoring instead a centrally concentrated cluster of stellar-mass black holes that accounts for 4.6% of the mass of ω Cen. Zocchi et al. (2017) and Zocchi et al. (2019) stress that the presence of an IMBH in ω Cen can be degenerate with orbital anisotropy and with the presence of a cluster of mass-segregated stellar-mass black holes. Evans et al. (2022) concluded that a core of stellar remnants can account for up to ~5 × 10⁵ M_⊙ of the inner dark mass derived in their analysis, with dark matter also being a potential explanation. Therefore, the structure and composition of the mass content of ω Cen remains a subject of ongoing debate.

In recent years, MSP timing observations have been used to probe the gravitational potential and mass contents of GCs, with a growing wealth of data being produced¹ and analyzed in this context (Prager et al. 2017; Freire et al. 2017; Perera et al. 2017a,b; Kızıltan et al. 2017a,b; Abbate et al. 2018, 2019; Xie & Wang 2020; Ridolfi et al. 2021; Wang & Xie 2021; Zhang et al. 2023; Dai et al. 2023; Padmanabh et al. 2024; Xie et al. 2024; Vleeschower et al. 2024). We also refer the interested reader to earlier studies where this methodology was explored, such as in Phinney (1993) and Anderson (1993).

In Dai et al. (2023), constraints on accelerations from five centrally located (with projected radii of within a few parsecs of the center) MSPs in ω Cen were measured via Doppler-induced effects in timing observations from the pulsars discovered in Dai et al. (2020). Dai et al. (2023) reported a discrepancy between the extremal bounds derived from the baryonic profile used to model the mass content of ω Cen and the comparatively high inferred pulsar accelerations. This suggests that MSP line-of- sight (LOS) accelerations could offer information to bolster our understanding of the structure of ω Cen, adding to the previous discussion. In particular, robust constraints on accelerations near the central region could allow us to discern between the various profiles that have been put forward with unprecedented accuracy, with the ability to identify sufficiently centrally concentrated distributions should they be present.

With the above in mind, in this work we carry out a Jeansbased analysis, exploiting velocity dispersion (proper motion and LOS) and surface brightness profile data with an implementation based on the publicly available code GravSphere² (Read & Steger 2017; Read et al. 2018; Genina et al. 2020; Collins et al. 2021). At the same time, we also include novel constraints on both the pulsar distribution and LOS acceleration data from Dai et al. (2023) as a self-consistent modification of the original likelihood function used in GravSphere. This allows us to fully exploit both the stellar and MSP timing-derived kinematics³.

An additional novel feature of our approach is that we decompose the inner mass distribution within ω Cen into a potential IMBH and an additional population tracing the observed distribution of MSPs. This allows us to break degeneracies between the mass in remnants and the possible presence of an IMBH, and provides a robust constraint on the spatial distribution of MSPs. This latter can also be used to shed light on the origin of the gamma-ray emission observed in ω Cen (Brown et al. 2019; Reynoso-Cordova et al. 2021; Dai et al. 2023, 2020).

This paper is organized as follows. In Sect. 2, we describe the mass modeling method that we use in this work and our treatment of the stellar kinematics. In Sect. 3, we discuss our stellar photometric and kinematic data. In Sect. 4, we discuss the pulsar timing data that we use to determine the pulsar accelerations, while in Sect. 5 we address their implementation in our methodology. In Sects. 6–8, we present our key results and we discuss them in the context of previous studies in the literature on ω Cen. Finally, in Sect. 9 we present our conclusions.

2 Jeans modeling of stellar kinematics

Our mass modeling method is adapted from the Jeans equations solver GravSphere. We start with a brief overview of how the latest version of GravSphere works. A more detailed discussion can be found in Read & Steger (2017), Genina et al. (2020), Collins et al. (2021), and Júlio et al. (2023).

In essence, GravSphere is designed to solve the Jeans equations under the assumption of a non-rotating, spherically symmetric stellar system which is in a steady state. This can be used to map velocity dispersion and surface brightness profile of collisionless tracer particles, such as stars in self-gravitating objects like galaxies and GCs, to the total underlying mass distribution.

We start by considering the dynamical equation governing the phase-space distribution function ƒ(x, v, t) of a collisionless system of particles, for example stars, given by the collisionless Boltzmann equation: $$\frac{\text{d}f}{\text{d}t}=\frac{\partial f}{\partial t}+{\nabla }_{x}f\cdot v-{\nabla }_{v}f\cdot {\nabla }_x\text{Φ}=0,$$(1)

where Φ is the gravitational potential resulting from all mass components in the system, and satisfying the Poisson equation: $${\nabla }_x^2\text{Φ}=4\pi G\rho ,$$(2)

where ρ is the total mass density of the system.

Although the distribution function ƒ is normally not directly accessible via observations, one can multiply Eq. (1) by powers of the velocity components and integrate in velocity space to arrive at a dynamic equation in terms of velocity moments which do relate to observable properties (e.g., see Binney & Tremaine 2008; Battaglia et al. 2013). In particular, taking the steady state solution $\frac{\partial f}{\partial t}=0$ for the first moment and multiplying by each velocity component v_j, one obtains a set of three well- known equations: the Jeans equations (Jeans 1922), which can be expressed as: $${\displaystyle \sum _{i=1}^3\frac{\partial {\nu }_{\star }\langle {\text{ν}}_i{\text{v}}_j\rangle }{\partial {\text{x}}_i}+v_{\star }\frac{\partial \text{Φ}}{\partial {\text{x}}_j}=0,}$$(3)

where ν_⋆(x) = ∫ d³vƒ is the stellar (tracer particle) number density and 〈v_iv_j〉 = ∫ d³vv_iv_jƒ correspond to the two-point velocity moments, with the brackets “〈〉” generally denoting an integral with the distribution function over velocity space.

For the case of a spherically symmetric distribution, one finds that the only non-trivial equation becomes that involving the radial component, which can be expressed in terms of the radial and tangential velocity dispersions (respectively): $${\sigma }_{r/t}^2=\langle {\text{ν}}_{r/t}^2\rangle -{\langle {\text{ν}}_{r/t}\rangle }^2,$$(4)

to give the spherical Jeans equation: $$\frac{1}{{\nu }_{\star }}\frac{\partial {\nu }_{\star }{\sigma }_r^2}{\partial r}+\frac{2\beta {\sigma }_r^2}{r}=-\frac{GM}{r^2},$$(5)

where $$\beta \equiv 1-\frac{{\sigma }_t^2}{{\sigma }_r^2}$$(6)

is the so-called velocity anisotropy⁴ and we adopt a Newtonian potential with enclosed mass M(r) within the 3D radial coordinate r. It should also be understood that ν_⋆,β, σ_r, σ_t are all, in general, functions of r.

Eq. (5) can be solved for the quantity ${\nu }_{\star }{\sigma }_r^2$ yielding the radial velocity dispersion (Mamon & Łokas 2005; van der Marel 1994): $${\sigma }_r^2(r)=\frac{1}{{\nu }_{\star }(r)𝑔(r)}{{{\displaystyle \int }}^{\text{}}}_r^{\infty }\frac{GM\left(r^{\prime }\right){\nu }_{\star }\left(r^{\prime }\right)}{r^{\prime 2}}𝑔\left(r^{\prime }\right)\text{d}{r}^{\prime },$$(7)

where $$𝑔(r)\equiv \exp \left(2{{\displaystyle \int }}^{\text{}}\frac{\beta (r)}{r}\text{d}r\right).$$(8)

The LOS projection of ${\sigma }_r^2$, which is a commonly used observable (e.g., Strigari et al. (2007)), is then given by (Binney & Mamon 1982) $${\sigma }_{\text{LOS}}^2(R)=\frac{2}{{\text{Σ}}_{\star }(R)}{{{\displaystyle \int }}^{\text{}}}_R^{\infty }\left(1-\frac{R^2}{r^2}\beta (r)\right)\frac{{\nu }_{\star }(r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r,$$(9)

where Σ_⋆(R) corresponds to the stellar surface density at a projected radius R, and one multiplies by the geometrical factors corresponding to the contributions of the radial and tangential components along the line of sight, taking a density-weighted average along the line of sight.

Having discussed some of the basic elements of a Jeans analysis, we now address an issue which arises due to the fact that (even if the photometric profiles of Σ_⋆,ν_⋆ are well constrained from surface brightness data) there is a clear degeneracy between the anisotropy β and the enclosed mass M(r) in Eqs. (7) and (9). This is known as the M − β degeneracy (or also, by extension, the ρ − β degeneracy, where ρ(r) is the density) and has been extensively discussed in the literature (Merrifield & Kent 1990; Wilkinson et al. 2002; Łokas & Mamon 2003; de Lorenzi et al. 2009; Read & Steger 2017). This problem is a manifestation of the fact that the Jeans equations alone are not a closed system, meaning that they cannot unambiguously constrain these two quantities without more information being provided.

Following, for example, Watkins et al. (2013); Van der Marel & Anderson (2010), in order to resolve the mass-anisotropy degeneracy, one can obtain expressions for the projected proper motions along the projected radius R and the tangential component (respectively), these are given by: $${\sigma }_{\text{PM},\text{R}}^2(R)=\frac{2}{{\text{Σ}}_{\star }(R)}{{{\displaystyle \int }}^{\text{}}}_R^{\infty }\left(1-\beta (r)+\frac{R^2}{r^2}\beta (r)\right)\frac{{\nu }_{\star }(r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r$$(10)

and $${\sigma }_{\text{PM},\text{t}}^2(R)=\frac{2}{{\text{Σ}}_{\star }(R)}{{{\displaystyle \int }}^{\text{}}}_R^{\infty }(1-\beta (r))\frac{{\nu }_{\star }(r){\sigma }_r^2(r)r}{\sqrt{r^2-R^2}}\text{d}r.$$(11)

Furthermore, following Merrifield & Kent (1990); Richardson & Fairbairn (2014), to help eliminate the ρ − β degeneracy, one can integrate with higher moments of velocity in the collisionless Boltzmann equation, Eq. (1), to obtain two independent observables known as virial shape parameters (VSPs), which are given by $$\text{\hspace{0.17em}VSP1\hspace{0.17em}}=\frac{2}{5}{{{\displaystyle \int }}^{\text{}}}_0^{\infty }GM{\nu }_{\star }(5-2\beta ){\sigma }_r^{2}r\text{d}r={{{\displaystyle \int }}^{\text{}}}_0^{\infty }{\text{Σ}}_{\star }\langle {\text{ν}}_{\text{LOS}}^4\rangle R\text{d}R$$(12)

and $$\text{VSP}2=\frac{4}{35}{{{\displaystyle \int }}^{\text{}}}_0^{\infty }GM{\nu }_{\star }(7-6\beta ){\sigma }_r^2{r}^3\text{d}r={{{\displaystyle \int }}^{\text{}}}_0^{\infty }{\text{Σ}}_{\star }\langle {\text{ν}}_{\text{LOS}}^4\rangle R^3\text{d}R,$$(13)

where one can obtain data from observations corresponding to the right-hand side of Eqs. (12) and (13) that can therefore constrain β.

For the photometric tracer density profile, in this paper, we introduce the generic αβγ profile, which has been used in a variety of contexts (Hernquist 1990; Zhao 1996; Di Cintio et al. 2014; Freundlich et al. 2020; Zentner et al. 2022) due to its ability to reproduce a diverse range of distributions analytically. This is given by a double power-law model given as follows: $${\nu }_{\star }(r)=\frac{{\rho }_c}{{\left(r/r_c\right)}^{\gamma }{\left(1+{\left(r/r_c\right)}^{\alpha }\right)}^{(\gamma -\beta )/\alpha }},$$(14)

where we have introduced the three exponent variables α, β, γ, and the scale radius and density r_c,ρ_c, respectively. In previous versions of GravSphere, a summation of mass components given by Plummer sphere profiles (Plummer 1911) was used. Plummer models are commonly used in this context and have been found to provide adequate fits to photometric profiles in many cases, while having a simple analytical form. However, we found the model we used was able to reproduce the observed surface brightness profile significantly better, particularly at larger projected radii where the profile becomes shallower than what Plummer-based models allow.

The projected tracer surface density, is then given by: $${\text{Σ}}_{\star }(R)={{{\displaystyle \int }}^{\text{}}}_{-\infty }^{\infty }{\nu }_{\star }(r(l))\text{d}l,$$(15)

where, for carrying out the integral, l corresponds to the coordinate along the LOS so that $r(l)=\sqrt{R^2+l^2}$. Unlike the original version of GravSphere, to aid computational efficiency, we pre-compute this projection before fitting the kinematic and acceleration data, thereby assuming that the uncertainty on the photometric light profile is small as compared to the uncertainty on the velocity dispersion profiles. We explicitly checked that this is a good approximation for ω Cen.

To obtain the corresponding (cumulative) mass profile M_⋆(r), one can simply integrate Eq. (14) over space and multiply by a normalization factor needed to match the total mass of the distribution M_⋆. This factor corresponds to the mass-to-light ratio (if using luminosities rather than number densities or arbitrary units for the surface density profile), which is assumed to be constant throughout the profile.

For the full mass modeling, besides the photometric mass profile M_⋆(r), which accounts for the dynamical mass of the photometric distribution, including stars and any remnants or objects traced by this distribution, we add a central mass, M_cen(r), which emulates a generic central cluster of remnants. This central mass is modeled as a single Plummer sphere component (Plummer 1911) given by $$M_{\text{cen}}(r)=M_{\text{cen}}\frac{r^3}{r_{\text{cen}}^3}{\left(1+\frac{r^2}{r_{\text{cen}}^2}\right)}^{-3/2},$$(16)

where in this case we introduce the corresponding total mass and scale lengths M_cen and r_cen , respectively. To model the presence of an IMBH, we also introduce a point mass M_BH.

Lastly, since the full likelihood implementation of MSP data requires a model of the MSP distribution (see Sect. 5.2), we also include a mass component that follows the MSP profile. This allows us to consider the presence of a distribution of stellar remnants traced by the MSPs that are more centrally concentrated than the stars. This profile is given by $$M_p(r)=\frac{4M_p}{3\sqrt{\pi }}\frac{\text{Γ}[(1-\alpha )/2]}{\text{Γ}[-(\alpha /2+1)]}\frac{r^3}{r_0^3}{}_2{F}_1\left[\frac{3}{2},\frac{1-\alpha }{2};\frac{5}{2};-\frac{r^2}{r_0^2}\right],$$(17)

where ₂F₁ is the Gaussian hypergeometric function and M_p is the total mass of the distribution, which is defined (convergent) for the density exponent parameter α < −2, while r₀ is the length scale parameter. This equation is obtained by integrating the density with the functional form shown in Eq. (30) over space and multiplying it by the normalization factor to obtain the mass distribution. The resulting total mass profile is thus given by: $$M(r)=M_{\star }(r)+M_{\text{cen\hspace{0.17em}}}(r)+M_{\text{BH}}+M_p(r).$$(18)

While dark remnant models of ω Cen favor two-component profiles of light and heavy remnants (see discussion in Sect. 8.2 and references therein), in general, different populations of objects are expected to show different degrees of segregation based on their dynamical histories and their intrinsic masses. Being more massive than main-sequence stars, but less than heavy remnants (such as stellar-mass black holes), MSPs could indeed trace a distribution of intermediate concentration, should the system undergo sufficient mass segregation (e.g., Phinney 1993). Here we stress the role of pulsars is to act as tracers, rather than direct contributors to the kinematics. This mechanism is allowable, for instance, if other remnants of similar mass which are likely more abundant, such as white dwarfs and neutron stars, would undergo a similar level of mass segregation and thus be traced by the pulsars, even if these have negligible kinematic contributions by themselves. On the other hand, if more limited mass segregation has taken place, then these other remnants could be traced by the photometric profile instead (showing little segregation as with the lighter main-sequence stars) and only the heavy remnants (e.g., stellar-mass black holes) could segregate (as predicted by some models and simulations, see Sect. 8.2). In this case, a more concentrated distribution of pulsars could be explained by their mechanism of formation rather than mass segregation (see Sect. 7).

There is also no a priori reason to exclude the coexistence of an IMBH with a central remnant distribution, which is an important consideration given that both have been invoked to explain ω Cen’s kinematics. Therefore, to fully explore all these scenarios, their potential kinematic relevance, and to explicitly consider any potential degeneracies between them, we have decided to work with the generic multi-component model from Eq. (18)⁵.

We use the anisotropy profile assumed in GravSphere following the generalized form from Baes & van Hese (2007), which is an extension of the Osipkov-Merritt anisotropy profile (Osipkov 1979; Merritt 1985), and is given by $$\beta (r)={\beta }_0+\left({\beta }_{\infty }-{\beta }_0\right)\frac{1}{1+{\left(\frac{r_t}{r}\right)}^{\eta }},$$(19)

where β₀ is the anisotropy at the center, β_∞ its limit as it asymptotically approaches infinity, r_t is the transition scale between these two regimes, and η is the exponent which modulates the steepness of the transition.

3 Stellar kinematic and photometric data

We combine the LOS stellar kinematic data from Reijns et al. (2006) with Sollima et al. (2009), adjusting their individual line of sight velocities to match the value in Pechetti et al. (2024), (v_LOS〉 = 232.5 km/s. This accounts for any systematic offset between the different datasets (the shift is of order 1 km/s, which is smaller than the uncertainty on the individual stellar velocities). We perform a quality cut on these data, retaining only those stars with velocity error smaller than 4 km/s. This yields 1644 LOS velocities. We augment these with the Noyola et al. (2010) data that are pre-binned for three different choices of center; more on this, below. We take HST proper motion (PM) data from Bellini et al. (2017), augmenting this with Gaia proper motion data from Vasiliev & Baumgardt (2021). As with the LOS data, we perform a similar quality cut on the PM data, retaining stars with PM uncertainties smaller than 2 km/s at the distance of ω Cen (assumed to be 5.2 kpc; see Sect. 4), this leaves a combined total of ~200 000 PMs which pass this criterion (including also additional filters based on data quality and membership probabilities). The PM data are perspective corrected as in van de Ven et al. (2006), assuming a mean LOS velocity and distance, as above. The surface brightness data are taken from the Bellini et al. (2017) catalog and augmented with Gaia photometry from de Boer et al. (2019). We explore three different choices of center for ω Cen taken from Anderson & Van der Marel (2010) (RA, Dec = 201.69683333, −47.47956944), Noyola et al. (2008) (RA, Dec = 201.69184583, −47.47911111) and the kinematic center from Noyola et al. (2010) (RA, Dec = 201.69630208, −47.47835389).

These centers are self-consistently computed for the photometric light profile, LOS and PM stellar kinematic data.

The LOS velocity and PM data are then binned using the binulator, as described in Collins et al. (2021). We use 100 logarithmically spaced bins, ensuring that no bin has less than 100 stars in it. For the LOS data, we augment the post-binned data with the pre-binned data from Noyola et al. (2010), using the consistent center reported in that work. binulator calculates the two virial shape parameters, VSP1 and VSP2 (see Sect. 8.1) and their uncertainties. As such, we use these also in our fits to assist with breaking the mass-anisotropy degeneracy (see Sect. 2). We explicitly checked whether our choices of binning parameters impact our results. Replacing our binned PM data with the binned data from Watkins et al. (2015a) yielded similar results for our mass components, with the main difference being a somewhat more extended distribution for the central mass component, which had a comparable total mass, showing no significant differences for the IMBH component. The stellar kinematic data and the best-fit photometric parameters used in this study have been made publicly available⁶.

4 Pulsar timing observations

Although pulsars are known to be remarkably regular sources, over sufficiently long timescales, it is in fact possible to detect significant variations in the observed period of these objects. This effect encodes information on the relative acceleration between the observer and pulsar reference frames, as well as the intrinsic spin-down of the pulsar (e.g., see Phinney 1993), and can be calculated through the following equation: $${\left(\frac{\dot{P}}{P}\right)}_{\text{obs}}=\frac{a_{\text{LOS}}+a_{\text{S}}+a_{\text{g}}}{c}+{\left(\frac{\dot{P}}{P}\right)}_{\text{int}},$$(20)

where c is the speed of light and P_obs/int, Ṗ_obs/int denote the period and its time derivative in the observer / pulsar rest frames (respectively), with the second term on the right-hand side corresponding to the intrinsic spin-down contribution of the pulsar. The first term on the right-hand side accounts for the various contributions due to the relative motions between the pulsar and the observer which induce an effective time-dependent Doppler shift in the observed period. We will now address the various components which constitute this term.

a_𝑔 is the relative difference in LOS accelerations due to the differential Galactic rotations of the GC and the Solar System. Following Dai et al. (2023); Freire et al. (2017) and Section 5.1.2 of Prager et al. (2017), we use Equation (5) from Nice & Taylor (1995) yielding $$a_g=-\cos (b)\left(\frac{{\text{Θ}}_0^2}{R_0}\right)\left(\cos (\ell )+\frac{\vartheta }{{\sin }^2(\ell )+{\vartheta }^2}\right){\text{ms}}^{-2},$$(21)

where ϑ ≡ (d/R₀)cos(b) − cos(ℓ). R₀ = 8.34 ± 0.16 kpc is the distance from the Sun to the Galactic center and Θ₀ = 240 ± 8 km/s is the Galactic rotational speed at that point, both of which are obtained from Reid et al. (2014). The distance from the cluster is fixed to d = 5.2 kpc. This value was originally chosen for consistency with the values reported in Bellini et al. (2014); Watkins et al. (2015a,b), and was subsequently validated throughout our analysis (see Sect. 8.1). We use the Galactic coordinates ℓ = 309.1°, b = 14.97° corresponding to the center of ω Cen located at RA = 13:26:47.24, Dec = −47:28:46.5 adopted in Dai et al. (2023)⁷.

For the propagation of errors, we used the 68% confidence limit (CL) interval about the median obtained after running one million iterations, where a set of random samples is generated for the input parameters Θ₀ , R₀ , under the assumption that they are independent and normally distributed with median values and errors as cited.

a_S accounts for the proper motion of the GC which induces an apparent LOS acceleration contribution – the Shklovskii effect – which is given by (Shklovskii 1970; Dai et al. 2023) $$a_S=3.78\times {10}^{-12}\left(\frac{d}{5.2\text{\hspace{0.17em}kpc}}\right){\left(\frac{{\mu }_T}{\text{mas\hspace{0.17em}}{\text{yr}}^{-1}}\right)}^2\text{m\hspace{0.17em}}{\text{s}}^{-2},$$(22)

where ${\mu }_T=\sqrt{{\mu }_{\delta }^2+{\mu }_{\alpha \ast }^2}$ is the proper motion, and we take µ_δ = −6.7445 ± 0.0019 mas yr⁻¹, μ_α* = −3.1925 ± 0.0022 mas yr⁻¹ (Helmi et al. 2018).

a_LOS corresponds to the LOS component of the acceleration caused by the gravitational field of the GC on the pulsar and, under the assumption of spherical symmetry, is given by $$a_{\text{LOS}}(r,l)=-G\left(\frac{l}{r}\right)\frac{M(r)}{r^2},$$(23)

where l in this case is the longitudinal component (i.e., along the line of sight) of the position of the pulsar from the center of the GC and r is the total distance, so that $r=\sqrt{l^2+R^2}$, where R is the projection orthogonal to the line of sight and is typically the one that can be measured directly. M(r) denotes the enclosed mass of the spherical density distribution, while the geometric factor l/r projects the LOS component of the total radial acceleration experienced by the pulsar. This is the component of the observed period derivative which traces the mass contents of the GC, and is therefore the reason why these timing observations are of interest for the present study.

While the intrinsic spin-down component is subject to significant uncertainties and is difficult to measure directly, it can, however, be approximated to a model of magnetic dipole emission with breaking index n = 3, yielding (Prager et al. 2017; Dai et al. 2023) $$a_{\text{int}}\equiv c{\left(\frac{\dot{P}}{P}\right)}_{\text{int}}\approx 7.96\times {10}^{-10}{\left(\frac{B}{2\times {10}^8\text{G}}\right)}^2{\left(\frac{2\text{ms}}{P_{\text{obs}}}\right)}^2\text{m\hspace{0.17em}}{\text{s}}^{-2},$$(24)

where B is the surface magnetic field strength of the pulsar, whose effective values can be constrained from observations where this term is dominant (see Sect. 5.3).

5 Likelihood definition

5.1 Implementation of stellar kinematics constraints

For the purposes of fitting the stellar kinematics observables discussed in Sect. 2, one can construct a general log-likelihood function for the proper motion and LOS velocity dispersions, and the virial shape parameters, taking the form: $$\ln ℒ(\theta )=-\frac{1}{2}{\displaystyle \sum _y{\chi }_y^2},$$(25)

where we have the general chi-square variable for an observable y defined as $${\chi }_y^2\equiv {\displaystyle \sum _i\frac{{\left[y_{i,\text{obs}}-y_i(\theta )\right]}^2}{\delta y_i^2}},$$(26)

and observables are taken to be normally distributed with median values {y_{i, obs}} and corresponding standard errors {δy_i}, and {y_i(θ)} correspond to the predicted variables as a function of the mass and (symmetrized) anisotropy model parameters θ.

5.2 Implementation of pulsar constraints

There are various possible ways to implement the pulsar timing constraints. However, there are some subtleties that should be considered when doing so. For example, a pointwise optimization for each of the pulsars over the nuisance parameter l can be quite numerically expensive if this has to be done for each step of a fitting routine. Also, errors or confidence intervals in acceleration constraints are not originally given in Dai et al. (2023) and can be difficult to quantify, particularly when including the intrinsic spin-down contribution.

To allow for an effective implementation and with the above discussion in mind, we include in the fit the quantity: $$a_{\text{obs}}\equiv c{\left(\frac{\dot{P}}{P}\right)}_{\text{obs}}-a_g-a_S=a_{\text{LOS}}+c{\left(\frac{\dot{P}}{P}\right)}_{\text{int}},$$(27)

where the right-hand side follows directly from Eq. (20) and corresponds to the theoretical prediction from the GC mass model used, including the contribution from the intrinsic spindown in Eq. (24). This observable is well constrained in each of its components. The errors can be propagated following our approach in Sect. 4: we simulate one million samples assuming that input variables are independent and normally distributed, with the specified errors and median values taken from the literature (see Sect. 4). We checked explicitly that the resulting distributions of a_obs were close to normal and obtained similar values to those predicted by common standard error propagation techniques (i.e., applying the central limit theorem to linear expansions of the functions). The resulting median values and errors on a_obs, as well as some of the data of the MSPs published in Dai et al. (2023), are shown in Table 1.

The term to be added to the log-likelihood function takes the form: $$\ln {ℒ}_{p,\mathrm{LOS}}\left({\theta }_M;\left\{l_i\right\},\left\{B_i\right\}\right)=-{\displaystyle \frac{1}{2}}{\displaystyle \sum _i^{N_p}\frac{{\left(a_{\text{obs},i}-a\left({\theta }_M;l_i,B_i\right)\right)}^2}{\delta a_i^2}},$$(28)

where a(θ_M ; l_i, B_i) is the predicted LOS acceleration including intrinsic spin-down (right-hand side of Eq. (27)) and θ_M are the set of parameters characterizing the mass model. The parameters {B_i} and {l_i} correspond to the (unknown) magnetic field and longitudinal components of each of the N_p pulsars (5 in our case) while {δa_i} are the propagated errors of Table 1.

To perform the fit, we marginalize over the parameters {B_i} and {l_i}. The l_i are allowed to vary over a wide range within the boundaries of the GC. The reason B is marginalized over is that it has not been measured for this sample of MSPs. However, a reasonable idea of its range can be obtained based on known measurements performed on similar MSP populations. This is possible by measuring period variations for MSPs with no GC potential contributions, meaning that, after correcting for differential Galactic rotation and the Shklovskii effect, these observations act as direct tracers of the intrinsic spin-down contributions. For binary systems that do not experience significant orbital variability, Doppler-induced changes can also be studied in orbital period derivatives, with the added benefit of being independent of intrinsic spin-down effects (e.g., Freire et al. 2017; Prager et al. 2017). To date, of the 18 known MSPs discovered in ω Cen, only the 5 used in this study have measured spin period derivatives and the only known binary from this set has not had its orbital period derivative measured (Dai et al. 2023). Therefore, we resort to existing studies of comparable MSP populations to constrain the intrinsic spin-down components. In particular, Prager et al. (2017) found using data from spin- dominated MSP populations in the ATNF catalog (Manchester et al. 2005) that the inferred distribution of B was well approximated to be log-normal with median log₁₀ B [G] = 8.47 and standard deviation of 0.33 dex, which we adopt as a prior. This allows for an ample range which captures some of the uncertainty in the intrinsic spin-down, while also being representative of the values expected for MSP populations.

a_obs can be regarded as a strict upper bound on the acceleration, in the sense that the contribution due to intrinsic spin-down will always be negative. This translates to a lower bound on the enclosed mass when the accelerations are negative. Implementing this as a prior during the fits would be a somewhat simpler approach and has the advantage of not depending on the uncertainties in the intrinsic spin-down. However, it is likely less constraining and would require the potentially numerically expensive point-wise optimization over the nuisance parameters {l_i} to determine whether the bound is satisfied.

To fully exploit the MSP data and help further constrain parameters, we adopt an approach similar to the one used by Prager et al. (2017) and Abbate et al. (2018). We add a positional component to the likelihood function, adopting the “generalized King model” profile which is typically used to model MSP populations. This has a surface density given by (Lugger et al. 1995): $$n(R)=n_0{\left[1+{\left(\frac{R}{r_0}\right)}^2\right]}^{\alpha /2},$$(29)

where we have two additional free parameters: the exponent α, which modulates the steepness of the distribution; and the scale radius r₀. (We also have the additional parameter n₀, which is the central surface density. However, this is not included in the fit as it is a global factor and the likelihood components that we shall introduce are independent of it.)

As seen in Grindlay et al. (2002), the 3D density has the form: $$n(r)=f\left(\alpha ,r_0,n_0\right){\left[1+{\left(\frac{r}{r_0}\right)}^2\right]}^{(\alpha -1)/2},$$(30)

where f (α, r₀, n₀) is a normalization factor that can be determined using the relation of the projected surface density (as in Eq. (15)) so that $$n(R)={\displaystyle {\int }_{-\infty }^{\infty }n}\left(\sqrt{R^2+l^2}\right)\text{d}l,$$(31)

and one can verify that this expression yields Eq. (29) with $$f\left(\alpha ,r_0,n_0\right)=\frac{\Gamma [(1-\alpha )/2]}{\sqrt{\pi }\Gamma (-\alpha /2)}\frac{n_0}{r_0},$$(32)

where Γ denotes the Euler gamma function.

To quantify the positional probability to be added to the likelihood function, we follow a similar approach to Appendix D of Anderson (1993) and Phinney (1993); Prager et al. (2017); Abbate et al. (2018). We express the total positional probability density (for an individual pulsar) as: $$P\left(l,R\mid \alpha ,r_0\right)=P\left(l\mid R,\alpha ,r_0\right)P\left(R\mid \alpha ,r_0\right),$$(33)

where $$\begin{array}{ll}\hfill P\left(l\mid R,\alpha ,r_0\right)\text{dl}=& \frac{n\left(\sqrt{R^2+l^2}\right)\text{d}l}{{\displaystyle {\int }_{-\infty }^{\infty }n}\left(\sqrt{R^2+l^2}\right)\text{d}{l}^{\prime }}\hfill \\ \hfill =& \frac{\Gamma [(1-\alpha )/2]}{\sqrt{\pi }\Gamma (-\alpha /2)}\frac{1}{r_0}\frac{{\left[1+\left(R^2+l^2\right)/r_0^2\right]}^{(\alpha -1)/2}}{{\left[1+{\left(R/r_0\right)}^2\right]}^{\alpha /2}}\text{\hspace{0.17em}d}l,\hfill \end{array}$$(34)

and⁸ $$\begin{array}{ll}P\left(R\mid \alpha ,r_0\right)\text{d}R& =\frac{Rn(R)\text{d}R}{{\displaystyle {\int }_0^{\infty }{R}^{\prime }}n\left(R^{\prime }\right)\text{d}{R}^{\prime }}\hfill \\ \hfill & =-(\alpha +2)\frac{R}{r_0^2}{\left[1+{\left(\frac{R}{r_0}\right)}^2\right]}^{\alpha /2}\text{\hspace{0.17em}d}R,\hfill \end{array}$$(35)

where in Eq. (34) we divide the infinitesimal amount of pulsars expected in a 3D ring of radius R, width dR, and length dl by those in the infinite cylinder of equal radius and width, recovering Eq. (3.7) of Phinney (1993). Similarly, for Eq. (35) we do the 2D analog, with the surface element of an annulus of width dR and radius R, integrating the whole disk with the projected surface density to obtain the relevant normalization factor in the denominator.

The complete positional probability for a given pulsar is thus given by: $$P\left(l,R\mid \alpha ,r_0\right)=-\frac{\alpha +2}{\sqrt{\pi }}\frac{\Gamma [(1-\alpha )/2]}{\Gamma [-\alpha /2]}\frac{R}{r_0^3}{\left[1+\left(l^2+R^2\right)/r_0^2\right]}^{(\alpha -1)/2},$$(36)

yielding the relevant (i.e., non-constant) contribution to the loglikelihood function for N_p pulsars: $$\begin{array}{l}\ln {ℒ}_{\text{pos}}\left(\left\{l_i\right\},\alpha ,r_0\right)=\frac{\alpha -1}{2}{\displaystyle \sum _i^{N_p}\ln }\left(1+\left(R_i^2+l_i^2\right)/r_0^2\right)\\ +N_p\left(\ln \frac{\Gamma [(1-\alpha )/2]}{\Gamma [-\alpha /2]}+\ln [-(\alpha +2)]-3\ln r_0\right).\end{array}$$(37)

As Chen et al. (2023) also reported the positions of 13 additional MSPs, but without the relevant timing data that trace kinematics, one can still use these positions to constrain the MSP distribution without the need to introduce any additional free parameters. For this purpose, we focus only on the projected radius component of the likelihood (as l becomes a redundant parameter). This gives the following log-likelihood (for ${N^{\prime }}_p$ pulsars): $$\begin{array}{ll}\ln {ℒ}_{\text{pos},\text{R}}\left(\alpha ,r_0\right)=& \frac{\alpha }{2}{\displaystyle \sum _i^{N_p^{\prime }}\ln }\left[1+{\left(\frac{R_i}{r_0}\right)}^2\right]\hfill \\ & +N_p^{\prime }\left(\ln [-(\alpha +2)]-2\ln r_0\right).\hfill \end{array}$$(38)

We note that in our treatment of the likelihood implementation of these MSP observables we have been careful to include normalization factors that depend on the parameters α and r₀. While these are not always explicitly shown in previous works (and can sometimes be ignored), they are necessary to include whenever they are parameter-dependent and not merely global factors of the likelihood or, correspondingly, additive constants in the log-likelihood function.

5.3 Full likelihood and priors

The full log-likelihood can now be summarized as $$\begin{array}{l}\ln {ℒ}_{\text{tot}}(\theta )=\ln {ℒ}_{\text{LOS}}+\ln {ℒ}_{\text{PM},\text{t}}+\ln {ℒ}_{\text{PM},\text{R}}+\ln {ℒ}_{\text{VSPI}}\\ +\ln {ℒ}_{\text{VSP}2}+\ln {ℒ}_{p,\text{LOS}}+\ln {ℒ}_{\text{pos}}+\ln {ℒ}_{\text{pos},\text{R}},\end{array}$$(39)

where θ includes all the model parameters which are simultaneously fitted, including anisotropy, photometric profile, MSP distribution, and mass model parameters. This is an important aspect of this study, as it allows for a fully consistent statistical treatment of the data and exploring potential degeneracies and correlations between different parameters, such as the aforementioned mass-anisotropy degeneracy. This interdependence of parameters may not otherwise be evident if one were to make separate fits or fix some of these parameters implicitly. Importantly, as will be shown in subsequent sections, it also allows for a full exploitation of the constraints imposed by the data in a way that differs from the results obtained from a separate treatment.

The prior on l is expressed in terms of the core-radius length scale r_c ≡ 4.54 pc of Baumgardt & Hilker (2018), which is also used by Dai et al. (2023). Following Read & Steger (2017), we also work with the symmetrized anisotropy, $$\tilde{\beta }\equiv \frac{\beta }{2-\beta },$$(40)

for $-1\le \tilde{\beta }\le 1$ to efficiently sample the (infinitely) broad range of possible values of β.

Table 2 shows the priors we use in our analysis, whose results will be discussed in the subsequent sections. The prior on α is based on the range of physically expected limits for pulsar populations, following the discussion in Phinney (1993); Prager et al. (2017)⁹. These limits are chosen to be broad and away from the tails of the posterior distribution in most cases, so that results are insensitive to them, while being sufficiently narrow to allow for satisfactory numerical efficiency. Priors in log-space are designed to efficiently span a broad range of several orders of magnitude.

The fits are performed using the nested sampling package dynesty (Speagle 2020; Koposov et al. 2022). This is based on nested sampling techniques (Skilling 2004, 2006) and the subsequently developed dynamic nested sampling technique (Higson et al. 2019), which we use, and is optimized for the inference of posterior distributions (as opposed to evidence estimation), using the bounding method from Feroz et al. (2009). The sampling method used is based on Neal (2003); Handley et al. (2015a,b). These methods are designed for more efficient sampling of complex posterior distributions which may present multimodalities and a potentially large number of dimensions in parameter space. The reason for this choice is that, in these settings, traditional Markov-chain Monte Carlo techniques tend to have less satisfactory performance, with numerical convergence being harder to achieve and the process being costly in CPU time. Another advantage of dynesty is that it is comparatively simple to implement, with the inclusion of an automated initialization routine and stopping criteria¹⁰.

6 Kinematic fits and constraints on a central mass

Figure 1 summarizes the fitting results obtained for the stellar kinematics and photometry observables used for the analysis. For the LOS and PM velocity dispersions and the two virial shape parameters, we obtain a best-fit reduced chi-square value of ${\chi }_{\nu }^2=1.61$ for our 4-parameter anisotropy profile and a 3- parameter mass model (stellar mass and Plummer mass and scale radius). We note that we have excluded both the black hole and MSP profile components due to the fact they do not have relevant kinematic contributions in our fits and adding them does not materially improve them¹¹. This is consistent with our finding that the data strongly favor an extended central mass in the inner regions and a stellar/photometric component in the outer ones, while excluding at 3σ CL an IMBH greater than 6 × 10³ M_⊙, or a kinematically relevant distribution of intermediate concentration traced by the MSPs. This bound on a putative IMBH mass is significantly more constraining than those reported in previous analyses. This result is illustrated in Fig. 2, which shows the posterior distributions of the mass model parameters and the mass and anisotropy profiles.

One reason for our bound being comparatively stringent is due to the extended central distribution that is favored, limiting the kinematic contribution of additional central components. This result was found to be statistically significant, since the ~10⁴ M_⊙ and ~10⁵ M_⊙ total mass values required for kinematically signficant IMBH and MSP profile contributions (respectively) exceed the 3σ limits of our derived posterior distributions, and are thus self-consistently excluded by the fit. We also found that, when fitting only the MSP component to model extended remnants, the MSP distribution was forced to emulate the concentration of our central mass favored by the kinematics, but failed to reproduce the more extended distribution of the MSPs, which is implicitly accounted for in the position likelihood components. This can be seen when comparing the cumulative distribution functions of the various components, which are examined in Sect. 7 (Fig. 4). This result indicates the need for a more concentrated mass component that is distinct from one traced by the pulsars.

Our results show a clear preference for a two-component model with an extended central mass of ~2–3 × 10⁵ M_⊙ and scale radius between 1.5 and 2.2 pc at the 3σ level. This is favored over an IMBH, leading to a 3σ upper bound of 6 × 10³ M_⊙ and placing a coexistence region for a putative IMBH of at most a few thousand solar masses. Our analysis also places a 3σ upper bound for the MSP distribution of 6 × 10⁴ M_⊙, strongly limiting the kinematic contribution of the more extended MSP component¹².

This is not present in previous analyses by construction, since these fit a single component of either a point mass or an extended distribution accounting for remnants, but not both simultaneously¹³. While previous studies have indicated the kinematic degeneracy between extended central distributions and an IMBH, in our analysis we are able to consider them simultaneously in a self-consistent fit.

One way the stellar kinematics favor this result is by having velocity dispersion profiles that are relatively flat in the central regions, with a value of ≲20 km/s within the inner ~0.1 pc, while still being sufficiently elevated within a few pc from the center to favor a significant extended central mass component of ~2–3 × 10⁵ M_⊙ that is needed in addition to the photometric one. Large (≳10⁴ M_⊙) IMBH models predict distinct cusps in velocity dispersion profiles (cf. Figure 5 of Zocchi et al. 2019), while more extended distributions do not. The fact that no such cusp is observed in our data, which has sufficient resolution to probe this central region, disfavors the presence of such component, and is instead consistent with the flatter profile predicted by the extended central mass component.

The similarity between the three velocity dispersion components explored is a manifestation of the approximate isotropy of the distribution, leading to a relatively well-constrained anisotropy profile. This can be observed clearly in Fig. C.2, where the three velocity dispersion components are plotted in the same figure, and is also apparent from the (symmetrized) anisotropy profile in Fig. 2.

The dynamical mass of the photometric component yields M_✶ = 2.68 ± 0.05 × 10⁶ M_⊙, which for a total luminosity of ~1.25 × 10⁶ L_⊙ would imply a mass-to-light ratio that is comparable, but arguably higher than those typical of stellar populations, showing consistency with previous studies, with some suggesting ω Cen to have high remnant fractions, for instance, in the form of white dwarfs (Dickson et al. 2024). For the photometric profile in Fig. 1 (middle right), we found this to be very well constrained by the data in morphology via the surface brightness data, yielding a very good fit with ${\chi }_{\nu }^2=1.13$ for the four-parameter model used. This, in addition to previous runs we performed where we fitted it with the kinematic data, justifies the simpler approach to fix the profile to its best-fit values.

Figure 3 illustrates the fits to the MSP kinematics corresponding to the LOS acceleration profiles for pulsars for the 5 MSPs with suitable timing solutions, noting that for the 13 remaining ones only projected positions were used as a constraint on the MSP distribution profile. We may now ask, more generally, the degree to which the MSP data are constraining of the overall mass distribution in this analysis. An important result in this regard is that our self-consistent fit with the MSP data yields a central mass component that is systematically more extended and massive when MSP accelerations are included, with in increase of ~20% in both the mass and scale radius of the distribution, with the black-hole and photometric components being comparatively unaffected (cf. Appendix B, where a fit is performed without including MSP accelerations).

The dominant contribution to this effect appears to be attributable to the most distant MSP acceleration data point, which effectively sets a lower limit on the enclosed mass at that point. The other data points were found to be independently consistent with the remaining pulsars, including the innermost one, whose relatively high acceleration was found to be potentially conflicting by Dai et al. (2023) with the models they explored, as these do not account for the presence of an inner dark mass such as the one favored in our analysis. Dai et al. (2023) also indicated a potential tension with this outermost point during their analysis. However, the degree to which such a tension is present in our analysis is necessarily dependent on the assumptions made on the intrinsic spin-down component, due to the fact that the lower bound implied by the observed apparent acceleration is within the allowable region of our analysis. While this data point favors a low intrinsic spin-down component, we found it to agree at the 2σ level with the observationally derived prior on the magnetic field strength parameter, with the other MSPs agreeing at the 1σ level. Considering also the significant observational and modeling uncertainties inherent to intrinsic spin-down determinations, we do not consider this to imply a significant tension in our analysis. It is interesting to note, however, that the constraining power effected by the outermost data point is robust with respect to intrinsic spin-down effects, as higher negative contributions to the acceleration would only increase the degree to which a more massive and extended distribution is favored.

In summary, we find the tension examined by Dai et al. (2023) to be alleviated by this analysis as a result of including a central dark mass component in our model, which emulates an extended cluster of stellar remnants (see discussion in Sect. 8.2), and by our statistical treatment of intrinsic spin-down modeling, which allows quantification of the significance of putative tensions.

For illustrative purposes, we also show in Fig. 3 the effect of concentrating ~15% of the central mass component in the form a 4 × 10⁴ M_⊙ IMBH, which is representative of some of the values in the literature (Noyola et al. 2008; Baumgardt 2017). We observe that, while in both models extremal LOS accelerations are allowable, the larger IMBH model would produce significant extremal accelerations over a larger region within the inner ~0.5 pc, comparable to the radius of the circle of influence r_inf ~ of the putative IMBH. A detection of such an extremal acceleration within this region would constitute a ‘smoking-gun’ signature of such an IMBH. Such scenario is disfavored by stellar kinematics in our analysis, except for the very inner region within ≲0.15 pc, where a smaller IMBH of at most a few thousand solar masses is still allowable¹⁴. This limited presence of extremal accelerations remains a falsifiable prediction in the context of our analysis, should new MSPs be detected in the innermost regions of ω Cen (cf. gray, dashed line of Fig. 3).

Fig. 1 Fits to stellar kinematic observables and the photometric light profile. The profiles are a function of projected radius at the assumed distance of 5.2 kpc. The green lines and bands correspond to the maximum posterior values, and the 68 and 95% CL centered at the median (respectively). For comparison, the gray (gray-dashed) line indicates the velocity dispersion obtained when a fraction of the dark component in the maximum posterior result is concentrated in the form of a 40 000 (10 000) M⊙ IMBH. Such IMBH models are strongly disfavored in our analysis. The stellar stellar kinematics datasets are described in Sect. 3. All data, including the photometric profile, have been self-consistently binned with the center from Anderson & Van der Marel (2010). Upper left: tangential proper motion velocity dispersion. Upper middle: radial proper motion velocity dispersion. Upper right: LOS velocity dispersion. Lower left: surface brightness profile used to determine the photometric component of the distribution. The red line shows the best-fit profile that was used throughout the analysis (see Sect. 2 for details). For presentation purposes, this profile has been normalized to match the central luminosity density of Trager et al. (1995) after correcting for extinction using the same approach as in Zocchi et al. (2017) with the reddening reported in Harris (1996) (2010 edition) catalog. Lower middle: posterior distribution for the virial shape parameter 1 as a result of the fits, with the red data point indicating the value with 1σ errors computed by binulator (see Sect. 3). Lower right: the same as the lower left figure, but for the virial shape parameter 2.

Fig. 2 Results for the mass - anisotropy profiles and mass model / pulsar distribution parameters. Upper left: 3D enclosed mass profile including photometric, central remnants, black hole, and MSP mass components. The solid lines correspond to the maximum posterior with 68 and 95% CL regions. The black and gray dashed lines indicate the upper limit of the 95% CL region for the black hole and MSP components, respectively. Upper right: symmetrized anisotropy profile including the maximum posterior and 68 and 95% CL regions. The posterior distributions for the anisotropy parameters are presented in Appendix A. The profile is close to isotropic, with a mild radial (tangential) anisotropy in the inner (outer) regions. This translates to the 1σ CL band being well within $\left|\tilde{\beta }\right|$, or |β| ≲ 0.15 for the vast majority of the range covered. Lower left: posterior distributions for the masses of the various fitted components of our analysis, indicating the median and 68% CL region at the top of each distribution. Lower right: posterior distributions for the morphological parameters of the fitted mass profiles with respective median and 68% CL regions.

Fig. 3 Absolute values of MSP LOS accelerations. The black solid line denotes the 95%CL upper bounds as inferred from the mass profiles of the posterior distribution. This region marks the boundary between the excluded region (red) and the allowable one (blue), where it is possible to find a value of l (i.e., a position within the cluster) compatible with the LOS acceleration at a given projected radius. The upwards (downwards) pointing triangles denote the observed MSP accelerations with positive (negative) values, as inferred from the MSP timing data. Errors are not included as they fall below the size of the data points. The vertical black lines denote the intrinsic spin-down components leading to the intrinsic LOS accelerations that trace the GC potential, based on the 68% CL regions in the magnetic field strength posteriors. This interval accounts for 84% of the posterior distribution as it includes the lower tail below the 16% percentile value in addition to the 68% CL interval contributions. The black, dotted line shows the corresponding bound for an illustrative model where ~15% of the central mass distribution is concentrated in the form of a 4 × 104 M⊙ IMBH. As can be seen, this model is still consistent with the pulsar accelerations (triangles), but is in tension with the proper motion and line of sight velocity data. A lower mass IMBH is still allowed (see the central allowable region that reaches to higher accelerations). The gray, dashed line indicates the projected position of the innermost detected MSP in ω Cen at ~0.86 pc from its center.

Fig. 4 Normalized cumulative distributions of MSPs (red), the more concentrated central mass emulating heavier stellar remnants (blue), stars (yellow), and X-ray-emitting sources (black and gray) in ω Cen as a function of projected radius. The inset shows a close-up view of the distribution normalized at a smaller radius, where all of the 18 MSPs are located. The inset includes the 68 and 95% CL regions of the MSP profile fit (shaded bands). The red dashed line denotes the median of the inferred cumulative distribution from the MSP 3D density from Eq. (30). The dashed-orange line indicates the predicted distribution derived from the stellar encounter rate that provides a remarkable match to the MSP distribution. We also show counts of X-ray sources observed in ω Cen studied by Henleywillis et al. (2018), showing the total count of 233 objects (black) and a subset of ~32 of the objects that share luminosities and X-ray colors compatible with known MSPs from other GCs (gray), as presented in Figure 10 of Heinke et al. (2005) (see also the discussion in Henleywillis et al. 2018). Over the radial range of the inset figure, this count is reduced to 105 and ~17 sources, respectively.

7 Analysis of the pulsar distribution

Figure 4 shows the (normalized) cumulative distributions of the MSP, photometric, and extended central mass components derived from our model fits. We also plot an independent dataset of X-ray sources shown for comparison purposes. The morphology of the MSP distribution is well constrained by the MSP data. It is clear from the figure that the MSP distribution is significantly more concentrated than the photometric component. The X-ray dataset was introduced to compare whether it may be independently traced by the components explored in this study. This is of particular interest given that 11 of the 18 known MSPs in ω Cen possess X-ray counterparts (Zhao & Heinke 2023; Dai et al. 2023). While it seems highly plausible that some of the MSPs are in fact found among these X-ray sources, interestingly, we find that the aggregate distribution of X-ray sources appears to be traced by the (more extended) stellar density profile, which, as implied by our previous discussion, was determined independently from stellar photometric data. This result appears to be essentially independent of which of the two X-ray sources is included, with the one sharing MSP-like properties (black) matching the total distribution (gray).

We now consider what the more concentrated MSP profile that we find here can tell us about MSP formation. MSPs have long been thought to descend from low mass X-ray binaries (e.g., see Bahramian et al. 2013; Verbunt & Freire 2014). These form dynamically in GCs at a rate set, to leading order, by the rate of stellar encounters¹⁵. Such models predict that the abundance of MSPs in GCs should scale with the (total) stellar encounter rate for the cluster Γ ∝ ʃ dr r²ρ_✶(r)²/σ(r), where ρ_✶(r) is the 3D density profile of stars, σ(r) is the 3D velocity dispersion profile of the stars, with ρ_✶(r)²/σ(r) being the encounter rate density, which encapsulates the rate of two-body stellar interactions (hence the ${\rho }_{*}^2$) and is enhanced by gravitational focusing (hence the inverse relation to σ). In practice, however, we may expect some deviation from this relation as the formation of MSPs depends on the binary stellar encounter rate that can differ from the single-star one (e.g., Verbunt & Freire 2014), binary population statistics, loss of stars from the cluster (impacted by the cluster escape velocity; Yin et al. 2024), neutron star retention fractions (e.g., Verbunt & Freire 2014), and the effects of mass segregation (e.g., see Bahramian et al. (2013) and references therein). Despite these caveats, numerous studies have reported a correlation between the MSP population abundances and the (total) stellar encounter rate in GCs (Bahramian et al. 2013; Zhao & Heinke 2022).

Previous studies have focused on establishing correlations between total encounter rates and abundances of MSPs in GCs. Following the above discussion, we decided to extend these analyses by assessing whether a similar relation exists at the intra-cluster level. In particular, we explore whether the MSP population density scales with the stellar encounter rate density so that ρ_p(r) ∝ ρ_✶(r)²/σ(r), which would naturally predict a linear correlation with the total encounter rate when both of these quantities are integrated over the volume of the cluster.

In the inset panel of Fig. 4, we compare our derived radial density profile of MSPs (shaded bands) with a model in which the MSP distribution depends solely on the stellar encounter rate (yellow dashed line). This gives an excellent match. Indeed, if we perform an indicative fit of the form ρ_p(r) ∝ ρ_✶(r)^γ/σ(r), we find that γ = 1.9 ± 0.3, showing excellent agreement for the central value of the inferred MSP distribution and well within its 68% CL bounds. Since we found that σ(r) ∝ ρ_✶(r)^0.2 to a good approximation for the region covered by the inset, this relation can also be expressed directly in terms of the stellar distribution as ρ_✶(r)^γ/σ(r) ∝ ρ_✶(r)^γ−0.2, so that ρ_p(r) ∝ ρ_✶(r)^1.7 for the central values.

Fig. 5 shows explicitly the dependence of the distributions considered as a function of the encounter rates, which for the case of ω Cen can be mapped to the profile via the enclosed encounter rate over a cylindrical volume with projected radius R, which we define as: $$\text{Γ}(<R)\propto {{{\displaystyle \int }}^{\text{}}}_{-\infty }^{\infty }\text{d}l{{{\displaystyle \int }}^{\text{}}}_0^R\text{d}{R}^{\prime }{R}^{\prime }{\rho }_{\star }{(r)}^2/\sigma (r),$$(41)

where $r=\sqrt{R^{\prime 2}+l^2}$.

Encounter rates are usually expressed in arbitrary units, which in our case have adapted to be consistent with Bahramian et al. (2013) for the total value (i.e., Γ(< ∞) ~ 90.4). For obtaining the constant of proportionality K mapping the encounter rate to the number of MSPs, we perform a least-squares fit (Fig. 5: red line) where we model the cumulative number of MSPs as: $$\text{\hspace{0.17em}Number\hspace{0.17em}of\hspace{0.17em}MSPs\hspace{0.17em}}(<R)=K\text{Γ}(<R)\text{,\hspace{0.17em}}$$(42)

yielding K ~ 0.25. This closely matches the central value of the total number of MSPs (~23) independently predicted by the phenomenological relation from Zhao & Heinke (2022) (Fig. 5: gray, dashed line). However, the 1σ upper bound of this relation implies more than double the abundance, suggesting that ω Cen could host a significant number of yet-undiscovered MSPs. One can also use the kinematic component of the MSP distribution to derive a strict upper bound on the MSP abundance assuming the totality of the dynamical mass being due to ~1.6 M_⊙ MSP binary systems, obtaining ~4 × 10⁴ at 3σ CL (c.f. Fig. 2).

The most significant result from Fig. 5 is the clear linear scaling observed for the MSP distribution, yielding a relation analogous to those from previous works for total encounter rates and abundances of GCs, but for the enclosed intra-cluster populations instead. This validates the model from Eq. (42), showing excellent consistency with the expectation from stellar encounter models of MSP formation discussed earlier.

On the other hand, for the case of the stellar density profile and the X-ray sources independently traced by it, we observe markedly different behavior, with a distinct nonlinear dependence followed in both cases. This further illustrates the differences between these distributions, with the linear scaling behavior being a unique characteristic of the MSP population.

To the best of our knowledge, while such scaling relations have been explored at the level of total encounter rates and MSP counts for GC populations, the present analysis is the first where a similar relationship is validated as a function of radius within a single GC, suggesting a functional form for the profile of the pulsar distribution within the cluster.

It is also interesting to note that, following the discussion of Sect. 8.2, since significant mass segregation for lighter stellar remnants (which MSPs are part of) is not expected for ω Cen, a different mechanism would be needed to explain their more concentrated distribution. This further validates the stellar encounter model, which naturally predicts the rate of encounters to be proportional to the square of the stellar density, and thus leads to a more concentrated distribution. This consideration, in conjunction with its abundance of recently discovered pulsars, makes ω Cen an ideal candidate to study such models of MSP formation¹⁶.

Fig. 5 MSP abundances as a function of stellar encounter rates. The red points are the observed cumulative number of MSPs as a function of the enclosed encounter rate at a projected radius R (upper x-axis). The red line corresponds to a least-squares linear fit, showing a clear linear dependence. The black points denote the same quantities for the full list of 233 X-ray sources from Henleywillis et al. (2018), while the orange line corresponds to the stellar distribution, normalized to match these sources. These show a distinct, nonlinear, dependence, with encounter rates that are not followed by the MSPs. The orange, dashed line is a linear extrapolation shown for comparison. Lastly, the gray dashed line indicates an upper estimate on the total number of MSPs as a function of total encounter rates of GCs, with 1σ error bands. This is based on the parametric fit performed by Zhao & Heinke (2022) using the phenomenological estimates of Bagchi et al. (2011), based on luminosity functions, and assuming the updated total stellar encounter rates of Bahramian et al. (2013). Units for ω Cen’s total encounter have been normalized to match the 90.4 central value of Bahramian et al. (2013).

8 Discussion

8.1 Method

The conflicting results in the literature regarding the nature of a central dark mass component in ω Cen are likely attributable to systematic differences in the data considered and the massanisotropy modeling adopted. For instance, Noyola et al. (2008) only considered a limited dataset of radial velocities for the stellar kinematics, while also assuming isotropy, which can lead to degenerate effects with anisotropic mass components (e.g., Zocchi et al. 2017; Read & Steger 2017). As noted in Sect. 8.2, degenerate effects with central remnants should also be accounted when considering the potential presence of IMBHs. These aspects have been extensively addressed in our analysis, with the inclusion of flexible mass-anisotropy modeling and a full self-consistent consideration multiple stellar kinematic data, including also MSP accelerations as an additional constraint.

Noyola et al. (2010) found a dependence of the inferred central mass with the assumed location of the kinematic center in their analysis, with the center used in Anderson & Van der Marel (2010) leading to a lower estimate of the inferred central mass. In a recent study from Pechetti et al. (2024), employing 3 independent center determinations, it was concluded that the favored center led to inferred stellar kinematics determinations that were consistent with the results from different datasets across multiple studies, including those of Anderson & Van der Marel (2010) using HST data. However, a significant tension was found with the originally assumed center in Noyola et al. (2008). This is also consistent with the findings of Anderson & Van der Marel (2010), where the Noyola et al. (2008) center was ruled out at high confidence and determined to be susceptible to systematic errors not originally considered in the study of Noyola et al. (2008).

We have addressed this issue using data with self-consistent centers (including our derived binning) for the kinematics and photometry, while also considering different centers such as the ones from Noyola et al. (2010) and Noyola et al. (2008). These results are shown in Appendix C (Figs. C.1 and C.2). While the inferred velocity dispersion profiles appear broadly consistent with those derived with the Anderson & Van der Marel (2010) center, we found that the central mass component derived with this center has a median mass (scale radius) that is ∼15% (∼11%) lower using the Noyola et al. (2010) center and ∼14% (∼10%) higher using the Noyola et al. (2008) center, with less significant effects on the remainder of the mass components in both cases. This indicates a somewhat more extended inferred central mass distribution when using the Noyola et al. (2010) center and a more concentrated one for the Noyola et al. (2008) center. This shows that the essential conclusions of our analysis are robust with respect to the choice of centers considered, with the main effect being a moderate variability in the concentration of the central mass component.

Another aspect in the literature concerns the various cluster distances adopted. We note that the close-to-isotropic distribution which is apparent in the data (cf. Figs. 2 and C.2) shows independent consistency with the 5.2 kpc adopted. This was further validated by performing a fit in which the distance was varied, showing a strong preference for the value adopted. This is not a trivial result given the different distance dependencies of PM and LOS velocity components and adds further support to the distance determinations in Watkins et al. (2015b), which exploit this effect. This value also shows agreement with the distance determinations of the Harris (1996) catalog (2010 edition) and the more recent ones from Baumgardt & Vasiliev (2021).

Van der Marel & Anderson (2010) reported central values for the central and outer anisotropy parameters of β₀ = 0.13 ± 0.02 and β_∞ = −0.53 ± 0.22, this is consistent at the 1-σ level with our determinations (cf. Fig. A.1), although our derived profile shows a higher degree of isotropy. This has a dependence on the cluster distance, which is significantly higher in our case (and in better agreement with recent determinations). Watkins et al. (2015a) found good agreement with isotropy for the central region of ω Cen, as measured by the proper motion anisotropy ratio, while finding some deviations for the outer regions of order ∼ 20%, indicating some radial anisotropy. Our proper motion data show comparable scatter, although the trend toward radial anisotropy at larger radii is less marked. In any case, the fact that we did not find strong differences when replacing our proper motion data with those of Watkins et al. (2015a) indicates that our mildly anisotropic models that are favored are also adequate in this case.

Another factor which may play a role in central mass determinations as seen, for instance, in Van der Marel & Anderson (2010), would be the photometric profile determinations. These authors found that different assumptions about the morphology of the photometric component led to differences in the degree to which central masses were favored. This is indicative of the importance of using a photometric profile which accurately reproduces the photometry. In this regard, the introduction of the αβγ profile as applied to this study, which was found to produce excellent fits to the photometry, constitutes a significant methodological improvement with respect to other modeling approaches which may nonetheless work well for other cases¹⁷.

To account for the presence of rotation, we followed section 5.4.3 of Fritz et al. (2016) by using full second velocity moments, which incorporate potential non-vanishing mean velocities acting as quadrature terms, rather than pure velocity dispersions.

As an indicative assessment of the potential effects of axisymmetry in our models, we also performed an analysis comparing with the results of Evans et al. (2022) where HST stellar kinematics from Watkins et al. (2015a) were fit using JAM axisymmetric models (Cappellari 2008), which also include a treatment of rotation. Using this dataset and adopting a Gaussian component for the dark mass (as in Eq. (D.1)) and a constant anisotropy as in Evans et al. (2022)¹⁸. After making these adjustments with the mass modeling and data used (which excludes the stellar kinematic and pulsar timing data in our main analysis), we were able to recover the central values for both the mass of the dark and photometric components of Evans et al. (2022) within ~15% accuracy when comparing with our central values. This suggests a limited effect which may be attributable to axisymmetry and rotation, and is conservative in that these data did not include rotation signatures the way we accounted for them in our main analysis, and it does not account also for potential differences due to the photometric profiles adopted which may not be intrinsic to these effects¹⁹. This is instructive to compare with D’Souza & Rix (2013), who concluded that flattened models with rotation led to mass estimates for the cluster that were 10% lower, with a 4.2% difference in the estimate being attributable to rotation. We also note the discussion of Van der Marel & Anderson (2010), where it is argued that the effects of flattening and rotation should be limited in the central region of interest. Thus, while we intend to explore axisymmetric models in future works, we expect these to have only a moderate quantitative effect for the purposes of this study.

We have also considered the effects of different mass parametrizations for our dark component, finding that the essential conclusions of our analysis remain unchanged. This is explored in Appendix D.

Finally, we also performed an additional analysis where we allow for a potential underestimation of errors in the LOS velocity dispersion profile. This is to account for effects such as increased shot noise, which, despite being explicitly considered, for instance, in Noyola et al. (2010), may in some cases be prone to being underestimated. This was included in the form of a free parameter that corresponds to a fractional error relative to the centrally observed value which was added in quadrature when computing the likelihood. This had little effect on the central values, also favoring a central dark component over an IMBH. The main difference in this case was that some degeneracy between the concentrated dark component and the more extended one from the pulsars was present, meaning that the data could not fully discriminate between these. While the inclusion of LOS velocities is important in constraining our mass models, this further attests to the robustness of our results, even if a significant underestimation of errors were to be present. An analysis with additional datasets not originally included in this work will also be performed in a follow-up study.

8.2 A dark cluster of remnants in ω Cen

It is apparent from our results in Sect. 6 that the two kinematically relevant components, namely, the extended central mass and photometric profiles, are very well constrained by the data. In this context, it is instructive to compare our results with the analyses of Zocchi et al. (2019); Baumgardt et al. (2019); Dickson et al. (2024), where the presence of a cluster of stellar- mass black holes is determined to be a viable alternative to the IMBH hypothesis to account for the observed velocity dispersions in ω Cen. Stellar-mass black holes also arise naturally in dynamical simulations, where they have been found to play an important role in the core dynamics of GC analogs, and are predicted to be hosted in significant numbers by many presentday GCs (Kremer et al. 2020). Indeed, a concentrated cluster of remnants appears to be the only viable interpretation for the dominant mass component favored for the central region that our analysis favors over any of the other mass components considered. A two-component model of light and heavier remnants as the kinematically relevant contributions is indeed favored by our model, in agreement with the theoretically motivated assumptions of these studies.

Dickson et al. (2024) performed a recent analysis of ω Cen among other GCs, in which limepy distribution function models (Gieles & Zocchi 2015) were used to obtain realistic distributions for an ample spectrum of species comprising stellar-mass black holes, white dwarfs, neutron stars and main sequence stars. These were subsequently fitted to velocity dispersion and stellar mass function data following the procedure described in Dickson et al. (2023). A key finding from this study is the presence of a substantial central mass component dominated by stellar-mass black holes, with total mass ${1.82}_{-0.06}^{+0.05}\times {10}^5{\text{M}}_{\odot }$, while the remaining lighter remnants follow the same profile to a good approximation. This favors a two-component mass model whereby the segregation of lighter objects is inhibited. This shows broad agreement with previous studies of ω Cen where this has been argued, owing to its young dynamical age, long two-body relaxation time, and black hole population abundance, which has been found to be anti-correlated with mass segregation (Dickson et al. 2024). We were able to explicitly compare our results with the projected mass profiles from Dickson et al. (2024). We found that for the region beyond ~1 pc, the profile from (Dickson et al. 2024) closely matched a Plummer sphere with scale radius ~1.5 pc²⁰, compared to our derived value of 1.88 ± 0.10 pc. While this implies a profile that is systematically more concentrated than the one we derive, it agrees at the ~1σ level with the ~20% less massive and extended profile we obtained without the inclusion of MSP accelerations. This result is shown in Appendix B (Figs. B.1 and B.2), where a fit without the inclusion of MSP accelerations is performed for comparison. Interestingly, this suggests the analysis of Dickson et al. (2024) may benefit from the inclusion of MSP accelerations following our implementation, leading to a profile that would likely be more consistent with the one we obtain. Notwithstanding systematic uncertainties due to modeling assumptions, which should be fully accounted for a more definitive comparison, the similarity of these profiles, and the fact that lighter remnants would not be able to undergo the degree of segregation implied by our derived central mass distribution, stellar-mass black holes indeed appear to be the favored interpretation for our derived central component.

8.3 An IMBH in ω Cen?

Häberle et al. (2024) recently derived a lower bound on the mass of an IMBH in ω Cen of 8200 M_⊙ based on the observation of centrally located fast-moving stars, which is in apparent tension with our 3σ bound of 6 × 10³ M_⊙. However, this bound is dependent on the assumed escape velocity for the cluster, which, for instance, is ~10% higher in our analysis due to the presence of an extended dark mass attributable to stellar remnants (see Sect. 8.2).

Pechetti et al. (2024), on the other hand, found no evidence of an IMBH in the form of fast-moving stars in the inner ~0.5 pc region predicted by N-body IMBH simulations from Baumgardt (2017). Baumgardt et al. (2019) arrived at a similar conclusion based on the lack of observed fast-moving stars in their analysis, concluding that a ~4 × 10⁴ M_⊙ IMBH, favored by other analyses (Noyola et al. 2008; Baumgardt 2017), is ruled out in favor of a cluster of remnants, which is consistent with our results.

IMBH formation simulations in clusters have also found it challenging to produce black holes exceeding ~500 M_⊙ due to the ejection from gravitational wave recoils in black hole merger interactions (Arca Sedda et al. 2023; Fujii et al. 2024). Interestingly, in recent simulations from Fujii et al. (2024), it was found to be possible to produce ≳10³ M_⊙ IMBHs in GC analogs via very-massive star formation processes. However, the authors conclude that this formation scenario disallows IMBHs of masses greater than 10⁴ M_⊙, even for the case of very massive GCs such as ω Cen. This poses a theoretical challenge for any claims of IMBH detections in GCs exceeding this limit.

Many works have suggested that ω Cen is the stripped remnant of a nucleated dwarf (Bekki & Freeman 2003; Wirth et al. 2020; Johnson et al. 2020; Gray et al. 2024). If so, it should have a dark matter halo and a correspondingly higher escape velocity. This could explain the fast-moving stars found in Häberle et al. (2024) without the need for such a massive IMBH²¹, or also the observation from Pechetti et al. (2024) of a fast-moving star with a high (0.99) cluster membership probability that appears too distant from the center (≳1.5 pc) to be accounted for by an IMBH. To assess the viability of this explanation, more work is needed, including comparison with simulations and data not originally included in this analysis. This is, however, beyond the scope of this work. We will test this idea quantitatively in forthcoming papers in which we will add a dark matter halo to our ω Cen models.

9 Conclusion

The results and main findings of this work can be summarized as follows:

• We performed a combined analysis of stellar kinematics and LOS accelerations of MSPs to probe the mass content of ω Cen. We explored the existence of a central dark mass component and considered competing interpretations of its nature: a cluster of stellar remnants or an IMBH. We also exploited the data to model and constrain the MSP distribution, using this to probe MSP formation models, and as an additional mass component in the model that traces stellar remnants of intermediate mass.

• There are two key results from our analysis that are not present in previous works: One is that a two-component model with a significantly extended central mass distribution (with a scale radius of between 1.5 and 2.2 pc at the 3σ level) is strongly favored over both an IMBH (or more concentrated distributions) and the more extended component traced by the MSPs. The second result, which is a consequence of the first, is that we arrive at stringent bounds on the kinematic contributions of the two additional components considered, establishing a coexistence region between the favored central mass distribution and a putative IMBH with an upper mass limit of 6 × 10³ M_⊙ 3σ CL. Together, these findings favor a dark cluster of remnants rather than an IMBH as the explanation for the central kinematics of ω Cen.

• While the MSP profile – with a 3σ CL upper bound of 6 × 10⁴ M_⊙ – is found not to contribute meaningfully to the dynamical mass distribution of ω Cen, we find that the MSP timing constraints indeed play a significant role in the kinematics, favoring a central distribution that is ~20% more massive and extended.

• We find that the MSP density profile goes as: $${\rho }_p(r)\propto {\rho }_{\star }{(r)}^{\gamma }/\sigma (r),$$(43)

  with γ = 1.9 ± 0.3, which is consistent with models in which MSPs originate from stellar encounters. It is a natural expectation from this scenario that γ ~ 2. We also find that the stars and X-ray sources in ω Cen follow a similar density profile and are more radially extended than the MSPs. A significant new finding of our analysis is that the encounter rate–MSP abundance correlation can be extended to the intra-cluster level, motivating a profile for the spatial distribution of these objects.

• In future work, we will explore the presence of a dark matter halo in our mass model, as one is predicted to be present if ω Cen is the dissolved remnant of an accreted, nucleated dwarf (e.g., Gray et al. 2024).

• Despite limited observations leading to full timing solutions for only 5 of the 18 recently discovered MSPs in ω Cen, and more potential discoveries still underway (Chen et al. 2023), our analysis demonstrates the potential and promise of combining velocity dispersion data with pulsar accelerations to probe the inner mass distribution of GCs.

• Future observations may allow us to decipher whether or not extremal accelerations occur for central MSPs, which would be a “smoking-gun” signature of an IMBH, but is a disfavored scenario for ω Cen in our analysis.

With the number of MSP discoveries in GCs having recently doubled over a five-year period amid the advent of radio surveys with unprecedented sensitivity, such as those involving the MeerKAT telescope and the Square Kilometre Array (SKA), it is foreseeable that the methodology we introduce will be increasingly relevant for the understanding of the kinematics of these systems.

Acknowledgements

We thank the anonymous referee for very useful comments and feedback that have improved this paper. We are grateful to Renuka Pechetti for kindly providing us the photometric data used in this work. We would like to thank Vincent Hénault-Brunet, Giuseppina Battaglia, José María Arroyo Polonio, Paul Beck, Addy J. Evans, Jorge Sánchez Almeida, Claire S. Ye, and P.C.C. Freire for instructive discussions and their interest in our work. We also thank Joshua S. Speagle for providing us valuable information about the use of dynesty. We are grateful to Jorge Terol Calvo, Jorge García Farieta, and Ángel de Vicente for the invaluable support provided in using the high-performance computing systems at Insituto de Astrofísica de Canarias (IAC). We are grateful to Elena Pinetti and the team at the Cosmic Physics Center at Fermilab, including Alex Drlica-Wagner and Albert Stebbins, where this work was presented and discussed. The authors wish to acknowledge the contribution of the IAC High-Performance Computing support team and hardware facilities to the results of this research. A.B., F.C. and J.M.C. acknowledge funding received from the European Union through the grant “UNDARK” of the Widening participation and spreading excellence programme (project number 101159929). A.B. and J.M.C. acknowledge support from the MICINN through the grant “DarkMaps” PID2022-142142NB-I00. J.I.R. would like to acknowledge support from STFC grants ST/Y002865/1 and ST/Y002857/1. This work has made use of the following software packages: dynesty (Speagle 2020; Koposov et al. 2022), GravSphere (Read & Steger 2017; Read et al. 2018; Genina et al. 2020; Collins et al. 2021), corner.py (Foreman-Mackey 2016), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlid (Hunter 2007), Jupyter Notebook (Kluyver et al. 2016).

Appendix A Posterior distributions of the anisotropy parameters

Fig. A.1 Posterior distributions of the anisotropy parameters. The parameter space exhibits multiple degeneracies, but due to the correlations inherent to these, the resultant anisotropy profile is well constrained (cf. Fig. 2). We note that this is not inconsistent with our results and that these degeneracies are the result of our parametric modeling and not intrinsic to the anisotropy profile.

Appendix B Effects of excluding MSP accelerations

Fig. B.1 Posterior distributions for the masses of the components considered during a fit performed without the inclusion of MSP LOS accelerations into the likelihood function. The median for the central mass from the full fit (cf. Fig. 2) is ~ 20% greater than the one found during this fit, with the other components showing only moderate differences.

Fig. B.2 Posterior distributions for the morphological parameters of the mass profiles without the inclusion of MSP LOS accelerations into the likelihood function. The median for the central mass length scale from the full fit (cf. Fig. 2) is ~ 18% greater than the one found during this fit. The morphology of the MSP distribution is still constrained due to the projected radii included in the positional component of the likelihood function, favoring a somewhat different posterior distribution when accelerations are excluded. This, however, has a negligible effect on the remainder of the mass components in the absence of accelerations to constrain their kinematics.

Appendix C Use of different centers

Fig. C.1 Upper-left: Posterior distributions for the mass components considered in our kinematic models using the Noy10 center. The median of the central mass component for the And center, is ~ 15% lower, consistent with a more extended central mass distribution for the Noy10 center. The stellar mass obtained is only marginally smaller (≲ 3%), with the other IMBH and MSP components remaining subdominant without significant statistically meaningful differences. Upper-right: Same as the upper-left plot, but for the Noy08 center. In this case, the central mass for our main analysis with the And center is ~ 14% higher than the result shown, consistent with a more concentrated distribution for the Noy08 center, although still with the 1σ CL regions overlapping each other. A small increase of ≲ 2% is observed with respect to the stellar mass value, but these central values lie within their respective 1σ CL regions. The IMBH and MSP mass components, once more, show approximate consistency with the previous results.† Lower-left: Posterior distributions of the morphological parameters for the mass models using the Noy10 center. The median of the scale radius using the And center is ~ 11% lower, indicating a somewhat more extended central mass distribution for the Noy10 center, with the 1σ CL regions barely overlapping each other. The MSP parameters agree at the 1σ level. Lower-right: Same as the middle-left plot, but for the Noy08 center. In this case, the median of the scale radius using the And center is ~ 10% greater, indicating a somewhat more concentrated central mass distribution, but still showing 1σ CL compatibility. As with the Noy10 center, the MSP parameters show 1σ level consistency with the And center results. † Marginal differences in the morphology of the posterior distributions of the IMBH and MSP masses are not necessarily statistically meaningful, owing to the fact that, due to being inherently less constrained, these may require additional independent runs to fully establish such differences (which we did perform for the case of the And center). This, however, is not necessary for the purposes of the complementary analysis presented in this appendix.

Fig. C.2 Combined velocity dispersion profiles for the three components using the Anderson & Van der Marel (2010) center (And, left) Noyola et al. (2010) (Noy10, middle) and Noyola et al. (2008) (Noy08, right) kinematic centers. The maximum-posterior velocity dispersion profile (red) from our main analysis using the And center is shown for reference in all plots, yielding fits of comparable quality for all the cases considered. Due to the close-to-isotropic behavior of our inferred distribution at the 5.2 kpc distance employed, the three components of our maximum posterior fits are almost identical, with the large majority of the data points showing close overlap with each other.

Appendix D Use of different dark components

Figs. D.1 and D.2 explore the influence of differences in the adopted mass models for the extended dark component. In the first case, we adopt the functional form of the generalized profile in Eq. (17) with unconstrained priors (i.e., not ascribing to it constraints from the pulsar distribution as in our main analysis). This has the advantage of having a great deal of flexibility, as the exponent parameter α reproduces distributions of varying concentration, including the Plummer and King model profiles, and the parameter space explored allows anything from something approximating a point mass to a distribution of greater extension than the photometric profile. For both of these fits, we have left the point mass component due to a putative IMBH unchanged. In this case, we find that the distribution of all the mass components are compatible at the 1σ level, with an IMBH still being disfavored as a significant contributor to the kinematics. We note that the scale radius (which is mathematically identical to r₀ in Eq. (17)) may appear to indicate a more extended distribution when compared to our earlier result with the Plummer model, but this is due to a degeneracy in the r_cen − α_cen plane which is fixed in the former case and has little physical impact on our results. This result further corroborates that a) the favored dark component in our analysis is robust with respect to more general profile parametrizations that have the ability to reproduce it and b) that an additional extended distribution traced by the pulsars is not favored by these observations (consistent with the expected level of mass segregation in ω Cen, both theoretically and observationally).

Fig. D.1 Posterior distributions for the mass model parameters using the generalized mass profile in Eq. (17) for the extended dark component.

For the second case we adopt the density of a spherical Gaussian profile as in Evans et al. (2022), which can be expressed as: $${\rho }_{\text{cen,\hspace{0.17em}},\text{g}}(r)=\frac{M_{\text{cen\hspace{0.17em}}}}{{\left(2\pi r_{𝑔}^2\right)}^{3/2}}\exp \left(-\frac{r^2}{2r_{\text{g}}^2}\right).$$(D.1)

In this case, we observe a dark component that is systematically more concentrated, with a mass that is ~ 20% lower in this case, while the photometric profile is somewhat more massive, but with the 1σ intervals of the posteriors still overlapping. This, however, does not change the essential conclusions of our analysis of these distributions, namely, the presence of an extended dark component of ~ 2 − 3 × 10⁵ M_⊙ and the absence of a kinematically relevant IMBH component, capped at a few thousand solar masses.

Fig. D.2 Posterior distributions for the mass model parameters using a Gaussian profile (Eq. (D.1)) for the extended dark component.

References

All Tables

All Figures

Fig. 1 Fits to stellar kinematic observables and the photometric light profile. The profiles are a function of projected radius at the assumed distance of 5.2 kpc. The green lines and bands correspond to the maximum posterior values, and the 68 and 95% CL centered at the median (respectively). For comparison, the gray (gray-dashed) line indicates the velocity dispersion obtained when a fraction of the dark component in the maximum posterior result is concentrated in the form of a 40 000 (10 000) M⊙ IMBH. Such IMBH models are strongly disfavored in our analysis. The stellar stellar kinematics datasets are described in Sect. 3. All data, including the photometric profile, have been self-consistently binned with the center from Anderson & Van der Marel (2010). Upper left: tangential proper motion velocity dispersion. Upper middle: radial proper motion velocity dispersion. Upper right: LOS velocity dispersion. Lower left: surface brightness profile used to determine the photometric component of the distribution. The red line shows the best-fit profile that was used throughout the analysis (see Sect. 2 for details). For presentation purposes, this profile has been normalized to match the central luminosity density of Trager et al. (1995) after correcting for extinction using the same approach as in Zocchi et al. (2017) with the reddening reported in Harris (1996) (2010 edition) catalog. Lower middle: posterior distribution for the virial shape parameter 1 as a result of the fits, with the red data point indicating the value with 1σ errors computed by binulator (see Sect. 3). Lower right: the same as the lower left figure, but for the virial shape parameter 2.

In the text

Fig. 2 Results for the mass - anisotropy profiles and mass model / pulsar distribution parameters. Upper left: 3D enclosed mass profile including photometric, central remnants, black hole, and MSP mass components. The solid lines correspond to the maximum posterior with 68 and 95% CL regions. The black and gray dashed lines indicate the upper limit of the 95% CL region for the black hole and MSP components, respectively. Upper right: symmetrized anisotropy profile including the maximum posterior and 68 and 95% CL regions. The posterior distributions for the anisotropy parameters are presented in Appendix A. The profile is close to isotropic, with a mild radial (tangential) anisotropy in the inner (outer) regions. This translates to the 1σ CL band being well within $\left|\tilde{\beta }\right|$, or |β| ≲ 0.15 for the vast majority of the range covered. Lower left: posterior distributions for the masses of the various fitted components of our analysis, indicating the median and 68% CL region at the top of each distribution. Lower right: posterior distributions for the morphological parameters of the fitted mass profiles with respective median and 68% CL regions.

In the text

Fig. 3 Absolute values of MSP LOS accelerations. The black solid line denotes the 95%CL upper bounds as inferred from the mass profiles of the posterior distribution. This region marks the boundary between the excluded region (red) and the allowable one (blue), where it is possible to find a value of l (i.e., a position within the cluster) compatible with the LOS acceleration at a given projected radius. The upwards (downwards) pointing triangles denote the observed MSP accelerations with positive (negative) values, as inferred from the MSP timing data. Errors are not included as they fall below the size of the data points. The vertical black lines denote the intrinsic spin-down components leading to the intrinsic LOS accelerations that trace the GC potential, based on the 68% CL regions in the magnetic field strength posteriors. This interval accounts for 84% of the posterior distribution as it includes the lower tail below the 16% percentile value in addition to the 68% CL interval contributions. The black, dotted line shows the corresponding bound for an illustrative model where ~15% of the central mass distribution is concentrated in the form of a 4 × 104 M⊙ IMBH. As can be seen, this model is still consistent with the pulsar accelerations (triangles), but is in tension with the proper motion and line of sight velocity data. A lower mass IMBH is still allowed (see the central allowable region that reaches to higher accelerations). The gray, dashed line indicates the projected position of the innermost detected MSP in ω Cen at ~0.86 pc from its center.

In the text

Fig. 4 Normalized cumulative distributions of MSPs (red), the more concentrated central mass emulating heavier stellar remnants (blue), stars (yellow), and X-ray-emitting sources (black and gray) in ω Cen as a function of projected radius. The inset shows a close-up view of the distribution normalized at a smaller radius, where all of the 18 MSPs are located. The inset includes the 68 and 95% CL regions of the MSP profile fit (shaded bands). The red dashed line denotes the median of the inferred cumulative distribution from the MSP 3D density from Eq. (30). The dashed-orange line indicates the predicted distribution derived from the stellar encounter rate that provides a remarkable match to the MSP distribution. We also show counts of X-ray sources observed in ω Cen studied by Henleywillis et al. (2018), showing the total count of 233 objects (black) and a subset of ~32 of the objects that share luminosities and X-ray colors compatible with known MSPs from other GCs (gray), as presented in Figure 10 of Heinke et al. (2005) (see also the discussion in Henleywillis et al. 2018). Over the radial range of the inset figure, this count is reduced to 105 and ~17 sources, respectively.

In the text

Fig. 5 MSP abundances as a function of stellar encounter rates. The red points are the observed cumulative number of MSPs as a function of the enclosed encounter rate at a projected radius R (upper x-axis). The red line corresponds to a least-squares linear fit, showing a clear linear dependence. The black points denote the same quantities for the full list of 233 X-ray sources from Henleywillis et al. (2018), while the orange line corresponds to the stellar distribution, normalized to match these sources. These show a distinct, nonlinear, dependence, with encounter rates that are not followed by the MSPs. The orange, dashed line is a linear extrapolation shown for comparison. Lastly, the gray dashed line indicates an upper estimate on the total number of MSPs as a function of total encounter rates of GCs, with 1σ error bands. This is based on the parametric fit performed by Zhao & Heinke (2022) using the phenomenological estimates of Bagchi et al. (2011), based on luminosity functions, and assuming the updated total stellar encounter rates of Bahramian et al. (2013). Units for ω Cen’s total encounter have been normalized to match the 90.4 central value of Bahramian et al. (2013).

In the text

	Fig. A.1 Posterior distributions of the anisotropy parameters. The parameter space exhibits multiple degeneracies, but due to the correlations inherent to these, the resultant anisotropy profile is well constrained (cf. Fig. 2). We note that this is not inconsistent with our results and that these degeneracies are the result of our parametric modeling and not intrinsic to the anisotropy profile.
In the text

	Fig. B.1 Posterior distributions for the masses of the components considered during a fit performed without the inclusion of MSP LOS accelerations into the likelihood function. The median for the central mass from the full fit (cf. Fig. 2) is ~ 20% greater than the one found during this fit, with the other components showing only moderate differences.
In the text

Fig. B.2 Posterior distributions for the morphological parameters of the mass profiles without the inclusion of MSP LOS accelerations into the likelihood function. The median for the central mass length scale from the full fit (cf. Fig. 2) is ~ 18% greater than the one found during this fit. The morphology of the MSP distribution is still constrained due to the projected radii included in the positional component of the likelihood function, favoring a somewhat different posterior distribution when accelerations are excluded. This, however, has a negligible effect on the remainder of the mass components in the absence of accelerations to constrain their kinematics.

In the text

Fig. C.1 Upper-left: Posterior distributions for the mass components considered in our kinematic models using the Noy10 center. The median of the central mass component for the And center, is ~ 15% lower, consistent with a more extended central mass distribution for the Noy10 center. The stellar mass obtained is only marginally smaller (≲ 3%), with the other IMBH and MSP components remaining subdominant without significant statistically meaningful differences. Upper-right: Same as the upper-left plot, but for the Noy08 center. In this case, the central mass for our main analysis with the And center is ~ 14% higher than the result shown, consistent with a more concentrated distribution for the Noy08 center, although still with the 1σ CL regions overlapping each other. A small increase of ≲ 2% is observed with respect to the stellar mass value, but these central values lie within their respective 1σ CL regions. The IMBH and MSP mass components, once more, show approximate consistency with the previous results.† Lower-left: Posterior distributions of the morphological parameters for the mass models using the Noy10 center. The median of the scale radius using the And center is ~ 11% lower, indicating a somewhat more extended central mass distribution for the Noy10 center, with the 1σ CL regions barely overlapping each other. The MSP parameters agree at the 1σ level. Lower-right: Same as the middle-left plot, but for the Noy08 center. In this case, the median of the scale radius using the And center is ~ 10% greater, indicating a somewhat more concentrated central mass distribution, but still showing 1σ CL compatibility. As with the Noy10 center, the MSP parameters show 1σ level consistency with the And center results. † Marginal differences in the morphology of the posterior distributions of the IMBH and MSP masses are not necessarily statistically meaningful, owing to the fact that, due to being inherently less constrained, these may require additional independent runs to fully establish such differences (which we did perform for the case of the And center). This, however, is not necessary for the purposes of the complementary analysis presented in this appendix.

In the text

Fig. C.2 Combined velocity dispersion profiles for the three components using the Anderson & Van der Marel (2010) center (And, left) Noyola et al. (2010) (Noy10, middle) and Noyola et al. (2008) (Noy08, right) kinematic centers. The maximum-posterior velocity dispersion profile (red) from our main analysis using the And center is shown for reference in all plots, yielding fits of comparable quality for all the cases considered. Due to the close-to-isotropic behavior of our inferred distribution at the 5.2 kpc distance employed, the three components of our maximum posterior fits are almost identical, with the large majority of the data points showing close overlap with each other.

In the text

	Fig. D.1 Posterior distributions for the mass model parameters using the generalized mass profile in Eq. (17) for the extended dark component.
In the text

	Fig. D.2 Posterior distributions for the mass model parameters using a Gaussian profile (Eq. (D.1)) for the extended dark component.
In the text