LAPTH-038 24, CERN-TH-2024-131 /

## New constraints on the central mass contents of Omega Centauri from combined stellar kinematics and pulsar timing

Andrés Bañares-Hernández 1 2 , ⋆ , Francesca Calore 3 , Jorge Martin Camalich 1 2 4 , , , and Justin I. Read 5

- 1 Instituto de Astrofísica de Canarias, La Laguna, Tenerife, E-38200, Spain
- 2 Departamento de Astrofísica, Universidad de La Laguna
- 3 CNRS, Laboratoire d'Annecy-le-Vieux de Physique Théorique, 74940, Annecy, France
- 4 CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland
- 5 Department of Physics, University of Surrey, Guildford, GU2 7XH, United Kingdom

November 28, 2024

## ABSTRACT

Aims. We perform a combined analysis of stellar kinematics with line-of-sight accelerations of millisecond pulsars (MSPs) to probe the mass contents of Omega Centauri ( ω Cen). Our mass model includes the stellar mass distribution, a more concentrated mass component linked to the observed MSP population, a generic cluster of stellar remnants (assumed to be more concentrated than the stars and MSPs), and an intermediate-mass black hole (IMBH), allowing us to determine which of these is most favored.

Methods. We mass-model ω Cen using the package GravSphere to solve the Jeans equations, including constraints from proper motions, line-of-sight velocities, the surface density profile of the stars, the spatial distribution of MSPs, and the recently measured line-of-sight accelerations of a subset of these MSPs, self-consistently modeling their intrinsic spin-down. We explore the impact of di ff erent assumed centers of ω Cen on our results and we infer the posterior distributions of the model parameters from the combined likelihood using the nested sampling package dynesty .

Results. Our analysis favors an extended central mass of ∼ 2 -3 × 10 5 M over an IMBH, setting a 3 ⊙ σ upper limit on the IMBH mass of 6 × 10 3 M . We find that pulsar timing observations are an important additional constraint, favoring a central mass distribution ⊙ that is ∼ 20% more massive and extended than for models that are constrained by the stellar kinematics alone. Finally, we find a 3 σ CL upper bound of 6 × 10 4 M on the total mass traced by the MSPs, with the density profile following ⊙ ρ p ( r ) ∝ ρ⋆ ( r ) γ /σ ( r ) , with γ = 1 9 . ± 0 3 . , where ρ⋆ ( r ) is the stellar mass density and σ ( r ) is the stellar velocity dispersion profile. This favors models in which MSPs form via stellar encounters, such as the leading paradigm whereby MSPs are the progeny of low-mass X-ray binaries.

Conclusions. Our analysis demonstrates how combining stellar kinematics with MSP accelerations produces new constraints on mass models, shedding light on the presence or absence of IMBHs at the centers of globular clusters. Further, we provide the first validation of its kind where MSP positions are linked to their place of formation in globular clusters, showing excellent agreement with the expectation of stellar encounter models of MSP formation. This sets a promising precedent amid the advent of the rapidly growing number of observations and discoveries currently being experienced in this field.

Key words. galaxy / globular clusters: individual: Omega Centauri, stars / pulsars: general, stars: black holes, stars: kinematics and dynamics, galaxies: kinematics and dynamics, galaxies / star clusters: general, cosmology: dark matter

## 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 be strictly speaking a GC, but rather the remnant core of a tidally disrupted, nucleated, dwarf galaxy (Bekki &amp; 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 GaiaEnceladus 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 have their star formation first quenched by reionization and then reignited 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

⋆ a.banareshernandez@gmail.com

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 dark matter self-annihilation 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 &amp; 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 4 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 &amp; Anderson (2010), via a Jeansbased analysis, derived a 3 σ upper bound on the IMBH mass of 1 8 . × 10 4 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 5 M of the inner dark ⊙ mass derived in their analysis, with a 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 1 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; Abbate et al. 2019; Xie &amp; Wang 2020; Ridolfi et al. 2021; Wang &amp; 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 within a few pc from 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 tension 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 opens the possibility for MSP lineof-sight (LOS) accelerations to play a role in understanding the structure of ω Cen adding to the previous discussion. In partic-, ular, robust constraints on accelerations near the central region could allow one to discern between the various profiles postulated with unprecedented accuracy, with the ability to identify su ffi ciently 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 line-of-sight) and surface brightness profile data with an implementation based on the publicly available code GravSphere 2 (Read &amp; 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 line-of-sight acceleration data from Dai et al. (2023) as a self-consistent modification to the original likelihood function used in GravSphere . This allows us to fully exploit both the stellar and MSP timingderived kinematics. 3

1 For an updated census of observed MSPs in GCs, we refer the interested reader to this convenient web resource by P.C.C. Freire: https://www3.mpifr-bonn.mpg.de/staff/pfreire/GCpsr.html 2 https://github.com/justinread/gravsphere

3 We note that during the preparation of this manuscript we were made aware of a recent parallel study that builds on work in Hénault-Brunet et al. (2020) to also combine stellar kinematics with pulsar accelerations to improve mass models (Smith et al. 2024). Their approach is complementary to ours. They also include both stellar kinematic constraints and pulsar accelerations in a similar joint likelihood function, but they

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, as well as giving us 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 Sec. 2, we describe the mass modeling method that we use in this work and the treatment of the stellar kinematics. In Sec. 3, we discuss the stellar photometric and kinematic data we use. In Sec. 4, we discuss the pulsar timing data that we use to determine the pulsar accelerations, while in Sec. 5 we address their implementation in our methodology. In Secs. 6 - 8, we present our key results and we discuss them in the context of previous studies in the literature of ω Cen. Finally, in Sec. 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 &amp; 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 f ( x v , , t ) of a collisionless system of particles, e.g. stars, given by the collisionless Boltzmann equation:

$$\underset { \substack { n \text{ partic-} \\ l \text{ region} } } { \text{ n \text{ partic-} \\ \text{ s \text{ postu-} } } \quad \frac { d f } { d t } = \frac { \partial f } { \partial t } + \nabla _ { \text{x} } f \cdot \mathbf v - \nabla _ { \text{v} } f \cdot \nabla _ { \text{x} } \Phi = 0,$$

where Φ is the gravitational potential resulting from all mass components in the system, and satisfying the Poisson equation:

$$\inf _ { \substack { a \, \text{Jeans} ^ { - } \\ \text{tion and} } } \nabla _ { x } ^ { 2 } \Phi = 4 \pi G \rho,$$

where ρ is the total mass density of the system.

Although the distribution function f 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 &amp; Tremaine (2008); Battaglia et al. (2013)). In particular, taking the steady state solution ∂ f ∂ t = 0 for the first moment and multiplying by each velocity component v , one obtains a set of three wellj known equations: The Jeans equations (Jeans 1922), which can be expressed as:

$$\underset { \text{$r=r$ made} } { \overset { 3 } { \text{$r=r$ made} } } \sum _ { i = 1 } ^ { 3 } \frac { \partial \nu _ { * } \left \langle \nu _ { i } \nu _ { j } \right \rangle } { \partial x _ { i } } + \nu _ { * } \frac { \partial \Phi } { \partial x _ { j } } = 0,$$

implement this within the distribution function modeling code limepy (Gieles &amp; Zocchi 2015). They then apply their method to the 47 Tuc and Terzan 5 GCs, whereas we focus here on ω Cen.

where ν⋆ ( x ) = R d 3 v f is the stellar (tracer particle) number density and ⟨ v v i j ⟩ = R d 3 vv v i j f 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 / tangential velocity dispersions (respectively):

$$\sigma _ { r / t } ^ { 2 } & = \langle \nu _ { r / t } ^ { 2 } \rangle - \langle \nu _ { r / t } \rangle ^ { 2 } & ( 4 ) \quad \text{and}$$

to give the spherical Jeans equation:

$$\frac { 1 } { \nu _ { * } } \frac { \partial \nu _ { * } \sigma _ { r } ^ { 2 } } { \partial r } + \frac { 2 \beta \sigma _ { r } ^ { 2 } } { r } = - \frac { G M } { r ^ { 2 } }, & & ( 5 ) \quad \text{Frt} \\ & \quad \text{one can}$$

where

$$\beta & \equiv 1 - \frac { \sigma _ { t } ^ { 2 } } { \sigma _ { r } ^ { 2 } } & ( 6 ) & \text{given by}$$

is the so-called velocity anisotropy 4 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 ν⋆σ 2 r yielding the radial velocity dispersion (Mamon &amp; Łokas 2005; van der Marel 1994):

$$\sigma _ { r } ^ { 2 } ( r ) = \frac { 1 } { \nu _ { * } ( r ) g ( r ) } \int _ { r } ^ { \infty } \frac { G M ( r ^ { \prime } ) \nu _ { * } ( r ^ { \prime } ) } { r ^ { \prime 2 } } g ( r ^ { \prime } ) d r ^ { \prime }, \quad \quad ( 7 ) \quad \text{where } \text{$o1$} \\ \text{for $t$}$$

where

$$g ( r ) & \equiv \exp \left ( 2 \int \frac { \beta ( r ) } { r } d r \right ). & ( 8 ) & \text{jetty of } \cdot \\ & & \text{itty to rep}$$

The LOS projection of σ 2 r , which is a commonly used observable (e.g. Strigari et al. (2007)), is then given by (Binney &amp; Mamon 1982)

$$\sigma _ { \text{LOS} } ^ { 2 } ( R ) = \frac { 2 } { \Sigma _ { * } ( R ) } \int _ { R } ^ { \infty } \left ( 1 - \frac { R ^ { 2 } } { r ^ { 2 } } \beta ( r ) \right ) \frac { \nu _ { * } ( r ) \sigma _ { r } ^ { 2 } ( r ) r } { \sqrt { r ^ { 2 } - R ^ { 2 } } } d r, \quad \ \ ( 9 ) \quad \text{where} \ w \quad \text{and the} \quad \$$

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,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 &amp; Kent 1990; Wilkinson et al. 2002; Łokas &amp; Mamon 2003; de Lorenzi et al. 2009; Read &amp; 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.

4 Note that here we use the 1D tangential velocity dispersion. Some authors use the 2D tangential velocity dispersion, in which case the second term in Eq. (6) should be divided by a factor of two.

Following, for example, Watkins et al. (2013); Van der Marel &amp; 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:

$$\text{involving} \text{ } \text{$\text{$\sigma_{PMM,R}^{2}}(R) = \frac{2}{\Sigma_{*}(R)}\int_{R}^{\infty} \left(1-\beta(r) + \frac{R^{2}}{r^{2}}\beta(r) \right) \frac{\nu_{*}(r) \sigma_{r}^{2}(r) r} { \sqrt { r ^ { 2 } - R ^ { 2 } } } d r \ ( 1 0 )$$

and

$$\theta ^ { 2 } _ { \text{PM}, t } ( R ) = \frac { 2 } { \Sigma _ { * } ( R ) } \int _ { R } ^ { \infty } \left ( 1 - \beta ( r ) \right ) \frac { \nu _ { * } ( r ) \sigma _ { r } ^ { 2 } ( r ) r } { \sqrt { r ^ { 2 } - R ^ { 2 } } } d r.$$

Furthermore, following Merrifield &amp; Kent (1990); Richardson &amp; 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

$$\underset { 1 \text{coordi-} } { \text{wtonian} } \ \ V S P { 1 } = \frac { 2 } { 5 } \int _ { 0 } ^ { \infty } G M \nu _ { * } ( 5 - 2 \beta ) \sigma _ { r } ^ { 2 } r d r = \int _ { 0 } ^ { \infty } \Sigma _ { * } \left \langle \nu _ { L O S } ^ { 4 } \right \rangle R d R \ ( 1 2 )$$

and

$$\underset { \text{$g$ the ra$} } { \text{$r$} } _ { \text{$e$ Marel} } \ VSP 2 = \frac { 4 } { 3 5 } \int _ { 0 } ^ { \infty } G M \nu _ { * } ( 7 - 6 \beta ) \sigma _ { r } ^ { 2 } r ^ { 3 } d r = \int _ { 0 } ^ { \infty } \Sigma _ { * } \langle v _ { L O S } ^ { 4 } \rangle \, R ^ { 3 } d R,$$

where one can obtain data from observations corresponding to the RHS 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 _ { * } ( r ) = \frac { \rho _ { c } } { ( r / r _ { c } ) ^ { \gamma } ( 1 + ( r / r _ { c } ) ^ { \alpha } ) ^ { ( \gamma - \beta ) / \alpha } },$$

where we have introduced the three exponent variables α, β, γ, and the scale radius and density rc , ρ 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:

$$\ t h e \rho - \beta \\ \ t h e s i v e l y \quad \Sigma _ { * } ( R ) = \int _ { - \infty } ^ { \infty } \nu _ { * } ( r ( l ) ) \, d l, & & ( 1 5 )$$

where, for carrying out the integral, l corresponds to the coordinate along the LOS so that r l ( ) = √ R 2 + l 2 . Unlike the original version of GravSphere , to aid computational e ffi ciency, we precompute 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 }, & ( 1 6 ) & \beta ( r ) = \beta _ { \text{(} }$$

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 Sec. 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 { 4 M _ { p } } { 3 \sqrt { \pi } } \frac { \Gamma \left [ ( 1 - \alpha ) / 2 \right ] } { \Gamma \left [ - ( \alpha / 2 + 1 ) \right ] } \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 ], \, \begin{pmatrix} \text{$\text{of sight}$} \\ \langle v _ { L O S } \rangle = \\ \text{two then} \\ \text{smaller} \\ \text{We perf} \\ \dots \colon \dots \end{pmatrix}$$

where 2 F 1 is the Gaussian hypergeometric function and Mp is the total mass of the distribution, which is defined (convergent) for the density exponent parameter α &lt; -2, while r 0 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 _ { * } ( r ) + M _ { \text{een} } ( r ) + M _ { \text{BH} } + M _ { p } ( r ). \quad \ \ ( 1 8 ) \quad \text{assume} \quad \\ \quad \text{of } \sim \ \ 2 0$$

While dark remnant models of ω Cen favor two-component profiles of light and heavy remnants (see discussion in Sec. 8.2 and references therein), in general, di ff erent populations of objects are expected to show di ff erent 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 su ffi cient 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 Sec. 8.2). In this case, a more concentrated distribution of pulsars could be explained by their mechanism of formation rather than mass segregation (see Sec. 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). 5

We use the anisotropy profile assumed in GravSphere following the generalized form from Baes &amp; van Hese (2007), which is an extension of the Osipkov-Merritt anisotropy profile (Osipkov 1979; Merritt 1985), and is given by

$$( 1 6 ) \quad \beta ( r ) = \beta _ { 0 } + ( \beta _ { \infty } - \beta _ { 0 } ) \frac { 1 } { 1 + \left ( \frac { r _ { r } } { r } \right ) ^ { \eta } },$$

where β 0 is the anisotropy at the center, β ∞ its limit as it asymptotically approaches infinity, rt 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 o . / ff set between the di ff erent 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 di ff erent 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 &amp; 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 Sec. 4), this leaves a combined total of ∼ 200 000 PMs which pass this criterion (including , also additional filters based on data quality and mebership 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 di ff erent choices of center for ω Cen taken from Anderson &amp; 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

5 We note also that our choice of mass modeling is very general in that it makes few assumptions about the underlying distribution by allowing for multiple di ff erent mass components which are varied in nature and morphology. We also performed an additional analysis with di ff erent dark component parametrizations (Appendix D), finding that our main results are not a ff ected by these.

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 Sec. 8.1) and their uncertainties. As such, we use these also in our fits to assist with breaking the mass-anisotropy degeneracy (see Sec. 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 di ff erence being a somewhat more extended distribution for the central mass component, which had a comparable total mass, showing no significant di ff erences for the IMBH component. The stellar kinematic data and the best-fit photometric parameters used in this study have been made publicly available. 6

## 4. Pulsar timing observations

Although pulsars are known to be remarkably regular sources, over su ffi ciently long time scales, it is in fact possible to detect significant variations in the observed period of these objects. This e ff ect 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} }, & ( 2 0 ) \quad \stackrel { \sim } { \text{is the pr} } } \quad \text{is the pr} \\ & \quad \text{mass of}$$

where c is the speed of light and P obs int / , ˙ P obs int / denote the period and its time derivative in the observer / pulsar rest frames (respectively), with the second term in the RHS corresponding to the intrinsic spin-down contribution of the pulsar. The first term in the RHS accounts for the various contributions due to the relative motions between the pulsar and the observer which induce an e ff ective time-dependent Doppler shift in the observed period. We will now address the various components which constitute this term.

ag is the relative di ff erence in LOS accelerations due to the di ff erential 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 &amp; Taylor (1995) yielding

$$a _ { g } = - \cos ( b ) \left ( \frac { \Theta _ { 0 } ^ { 2 } } { R _ { 0 } } \right ) \left ( \cos ( \ell ) + \frac { \vartheta } { \sin ^ { 2 } ( \ell ) + \vartheta ^ { 2 } } \right ) \ m s ^ { - 2 }, \quad \ \ ( 2 1 ) \quad \text{where } \text{th}$$

where ϑ ≡ ( d R / 0 ) cos( b ) -cos( ℓ . ) R 0 = 8 34 . ± 0 16 kpc is the . distance from the Sun to the Galactic center and Θ 0 = 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 cho-. sen for consistency with the values reported in Bellini et al. (2014); Watkins et al. (2015a,b), and was subsequently validated throughout our analysis (see Sec. 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). 7

6 https://github.com/dadams42/OCenKinematics

7 This center was originally reported in Harris (1996) catalog (2010 edition). It is essentially equivalent to the Anderson &amp; Van der Marel (2010) center, which is the main center we adopt and present in this

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 Θ 0 , R 0 , under the assumption that they are independent and normally distributed with median values and errors as cited.

aS 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)

$$\overset { \prime \ y n u c } { \text{in differ} } _ { \text{e central} } \ a _ { S } = 3. 7 8 \times 1 0 ^ { - 1 2 } \left ( \frac { d } { 5. 2 \ k p c } \right ) \left ( \frac { \mu _ { T } } { \text{mas yr} ^ { - 1 } } \right ) ^ { 2 } m s ^ { - 2 }, \quad \ \ ( 2 2 )$$

where µ T = q µ δ 2 + µα 2 ∗ is the proper motion, and we take µδ = -6 7445 . ± 0 0019 mas yr . -1 , µα ∗ = -3 1925 . ± 0 0022 mas yr . -1 (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

$$\underset { \text{ation be} ^ { - } } { \text{le to de} ^ { - } } \\ \underset { \text{ation be} ^ { - } } { \text{le to de} ^ { - } }$$

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 = √ 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 di ffi cult 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)

$$\underset { \text{$\text{$System.} } { \text{$\text{$\text{$Sstm.} } } } { \text{$\text{$\text{$\text{$Sstm.} } } } } \ a _ { \text{int} } \equiv c \left ( \frac { \dot { P } } { P } \right ) _ { \text{int} } \approx 7. 9 6 \times 1 0 ^ { - 1 0 } \left ( \frac { B } { 2 \times 1 0 ^ { 8 } \, G } \right ) ^ { 2 } \left ( \frac { 2 m s } { P _ { \text{obs} } } \right ) ^ { 2 } \text{m} s ^ { - 2 }, \ \ ( 2 4 )$$

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

## 5. Likelihood definition

## 5.1. Implementation of stellar kinematics constraints

For the purposes of fitting the stellar kinematics observables discussed in Sec. 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:

$$\text{center of} \quad \text{in} \ \mathcal { L } ( \theta ) = - \frac { 1 } { 2 } \sum _ { y } \chi _ { y } ^ { 2 },$$

analysis for the stellar kinematic and photometric data. For the remaining centers (Appendix C), given the relatively extended positions of the MSPs and that these are not dependent on binning prescriptions, the choice of center should have a limited e ff ect on the results obtained.

where we have the general chi-square variable for an observable y defined as

$$\chi _ { y } ^ { 2 } & \equiv \sum _ { i } \frac { \left [ y _ { i, \, \text{obs} } - y _ { i } ( \Theta ) \right ] ^ { 2 } } { \delta y _ { i } ^ { 2 } }, & \text{trans} \quad & \text{tions. Fo} \\ & \text{bital vari} \\ & \text{period} \quad \text{period}$$

and observables are taken to be normally distributed with median values { yi , obs } and corresponding standard errors { δ yi } , and { yi GLYPH&lt;0&gt; θ GLYPH&lt;1&gt; } 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 di ffi cult to quantify, particularly when including the intrinsic spin-down contribution.

To allow for an e ff ective implementation and with the above discussion in mind, we include in the fit the quantity:

$$a _ { o b s } \equiv \, c \left ( \frac { \dot { P } } { P } \right ) _ { o b s } - a _ { g } - a _ { S } = a _ { L O S } + c \left ( \frac { \dot { P } } { P } \right ) _ { \text{int} }, \quad \text{(27)} \quad \text{tainties i} \\ \quad \text{training} \\ \quad \text{$\dots$} \quad \dots$}.$$

where the RHS follows directly from Eq. (20) and corresponds to the theoretical prediction from the GC mass model used, including the contribution from the intrinsic spin-down in Eq. (24). This observable is well constrained in each of its components. The errors can be propagated following our approach in Sec. 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 Sec.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 \mathcal { L } _ { p, \text{ LOS} } ( \boldsymbol \theta _ { M } ; \, \{ l _ { i } \}, \{ B _ { i } \} ) = - \frac { 1 } { 2 } \sum _ { i } ^ { N _ { p } } \frac { ( a _ { \text{obs}, i } - a ( \boldsymbol \theta _ { M } ; l _ { i }, B _ { i } ) ) ^ { 2 } } { \delta a _ { i } ^ { 2 } }, \ \ ( 2 8 ) \quad \text{introau} \quad \text{As $s$} \quad \text{for} \quad \text{$\lambda$ $s$} \quad \text{$\quad}$$

where a ( θ M ; l i , Bi ) is the predicted LOS acceleration including intrinsic spin-down (RHS of Eq. (27)) and θ M are the set of parameters characterizing the mass model. The parameters { Bi } and { l i } correspond to the (unknown) magnetic field and longitudinal components of each of the Np pulsars (5 in our case) while { δ ai } are the propagated errors of Table 1.

To perform the fit, we marginalize over the parameters { Bi } 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 di ff erential Galactic rotation and the Shklovskii e ff ect, 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 e ff ects (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 10 B [G] = 8 47 and standard devi-. ation 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 pointwise 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):

$$\text{obs where} \\ \text{redicted} \text{ applying } \ n ( R ) = n _ { 0 } \left [ 1 + \left ( \frac { R } { r _ { 0 } } \right ) ^ { 2 } \right ] ^ { \alpha / 2 },$$

where we have two additional free parameters: the exponent α, which modulates the steepness of the distribution; and the scale radius r 0 . (We also have the additional parameter n 0 , 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:

$$\text{including } \ n ( r ) = f ( \alpha, \ r _ { 0 }, \ n _ { 0 } ) \left [ 1 + \left ( \frac { r } { r _ { 0 } } \right ) ^ { 2 } \right ] ^ { ( \alpha - 1 ) / 2 },$$

where f ( α, r 0 , n 0 ) is a normalization factor that can be determined using the relation of the projected surface density (as in Eq. (15)) so that

$$\overset { \text{ters} } { \text{ithin the} } \ \ n ( R ) = \int _ { - \infty } ^ { \infty } n \left ( \sqrt { R ^ { 2 } + l ^ { 2 } } \right ) d l,$$

and one can verify that this expression yields Eq. (29) with

$$\begin{smallmatrix} \i \text{ known} \\ \i \text{ in GC} \end{smallmatrix} \begin{smallmatrix} f ( \alpha, \, r _ { 0 }, \, n _ { 0 } ) = \frac { \Gamma [ ( 1 - \alpha ) / 2 ] } { \sqrt { \pi } \Gamma ( - \alpha / 2 ) } \frac { n _ { 0 } } { r _ { 0 } }, \end{smallmatrix}$$

Table 1. Data from MSPs published in Dai et al. (2023) with the inclusion of our derived median values and 68% CL errors of a obs , as defined in Eq. (27), as well as the central values of R (see discussion following Eq. (23) for details).

| Name        | P obs [ms]         | ˙ P obs [10 - 20 s s - 1 ]   | RA [J2000]        | DEC [J2000]      |   R [pc] | a obs [10 - 9 ms - 2 ]   |
|-------------|--------------------|------------------------------|-------------------|------------------|----------|--------------------------|
| J1326-4728A | 4.108786192190(1)  | 2.738(2)                     | 13:26:39.6699(2)  | -47:30:11.641(3) |     2.89 | 1 . 914 ± 0 . 011        |
| J1326-4728B | 4.791869161014(2)  | -5.433(4)                    | 13:26:49.5688(3)  | -47:29:24.889(4) |     1.14 | - 3 . 483 ± 0 . 011      |
| J1326-4728C | 6.867859692327(6)  | 0.98(1)                      | 13:26:55.2219(6)  | -47:30:11.753(9) |     2.96 | 0 . 344 ± 0 . 012        |
| J1326-4728D | 4.578833468410(2)  | -4.110(4)                    | 13:26:32.7130(2)  | -47:28:40.053(3) |     3.72 | - 2 . 775 ± 0 . 012      |
| J1326-4728E | 4.2077170407405(4) | 1.628(1)                     | 13:26:42.67844(7) | -47:27:23.999(1) |     2.38 | 1 . 076 ± 0 . 011        |

Notes. Error propagation follows the approach outlined in Sec. 4.

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 ( l, R | \alpha, r _ { 0 } ) = P ( l | R, \alpha, r _ { 0 } ) P ( R | \alpha, r _ { 0 } ), & & ( 3 3 ) & & +$$

where:

$$\text{where} \colon & \\ P ( l | R, \alpha, r _ { 0 } ) \, d l = \frac { n \left ( \sqrt { R ^ { 2 } + l ^ { 2 } } \right ) d l } { \int _ { - \infty } ^ { \infty } n \left ( \sqrt { R ^ { 2 } + l ^ { 2 } } \right ) d l ^ { \prime } } & & \text{Since Cl} \\ = \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 + ( R / r _ { 0 } ) ^ { 2 } \right ] ^ { \alpha / 2 } } \, d l, \ ( 3 4 ) } & \\ & \ln \, \mathcal { C } \quad.$$

$$= \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 + ( R / r _ { 0 } ) ^ { 2 } \right ] ^ { \alpha / 2 } } \, d l, \ \ ( 3 4 )$$

and: 8

$$P ( R | \alpha, r _ { 0 } ) \, d R & = \frac { R n ( R ) \, d R } { \int _ { 0 } ^ { \infty } R ^ { \prime } n ( R ^ { \prime } ) \, d R ^ { \prime } } \quad & \text{Note tha} \\ & = - ( \alpha + 2 ) \frac { R } { r _ { 0 } ^ { 2 } } \left [ 1 + \left ( \frac { R } { r _ { 0 } } \right ) ^ { 2 } \right ] ^ { \alpha / 2 } d R, \quad & ( 3 5 ) \quad & \text{ization f} \\ \text{where in $E_{\alpha}$} \, / 2 \lambda \, \text{via indica the inpatient amount of matrices} \quad \text{over the} }$$

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 ( l, R | \alpha, r _ { 0 } ) = - \frac { \alpha + 2 } { \sqrt { \pi } } \frac { \Gamma \left [ ( 1 - \alpha ) / 2 \right ] } { \Gamma \left [ - \alpha / 2 \right ] } \frac { R } { r _ { 0 } ^ { 3 } } \left [ 1 + ( l ^ { 2 } + R ^ { 2 } ) / r _ { 0 } ^ { 2 } \right ] ^ { ( \alpha - 1 ) / 2 }, \quad \text{where $\theta$ mostly} \\ \text{distributi} \text{ aspect o!} }$$

8 Note that, in our definition, R is non-negative, so that we only integrate from 0 outwards (extending this to negative values would only change the integral by a global factor of 2, which would not a ff ect our results).

yielding the relevant (i.e., non-constant) contribution to the loglikelihood function for Np pulsars:

$$\text{media } D \\ \text{media } ( \text{2017} ) ; \\ \text{obability} \quad \ln \mathcal { L } _ { \text{pos} } ( \{ l _ { i } \}, \alpha, r _ { 0 } ) = \frac { \alpha - 1 } { 2 } \sum _ { i } ^ { N _ { p } } \ln ( 1 + ( R _ { i } ^ { 2 } + l _ { i } ^ { 2 } ) / r _ { 0 } ^ { 2 } ) \\ ( \text{33} ) \quad \quad + N _ { p } \left ( \ln \frac { \Gamma \left [ ( 1 - \alpha ) / 2 \right ] } { \Gamma \left [ - \alpha / 2 \right ] } + \ln [ - ( \alpha + 2 ) ] - 3 \ln r _ { 0 } \right ). \quad ( \text{37} ) \\ \text{c} \cdots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \cdots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \cdots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \, \mathcal { m } \, \dots \</dd>$$

Since Chen et al. (2023) have 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 Np ′ pulsars):

$$d l, \ ( 3 4 ) & \quad \cdot _ { \substack { \eta \\ \ln \mathcal { L } _ { \text{pos}, R } ( \alpha, r _ { 0 } ) } } ^ { \mu } & = \frac { \alpha } { 2 } \sum _ { i } ^ { N ^ { \prime } _ { p } } \ln \left [ 1 + \left ( \frac { R _ { i } } { r _ { 0 } } \right ) ^ { 2 } \right ] \\ & \quad + N ^ { \prime } _ { p } ( \ln [ - ( \alpha + 2 ) ] - 2 \ln r _ { 0 } ). \\ & \quad \cdot _ { \substack { \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \nu \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \colon \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \Pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \p \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pn \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \nop \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \right \nop \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi\pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi\pi \nu \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi\pi \p \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi\pi \nop \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \p \pi \nu \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \p \pi \p \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \p \pi \n \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \r \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \text \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pri \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \pi \type \pi \pi \pi \pi \pi$$

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 0 . 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

$$\text{ar is thus} \quad \ln \mathcal { L } _ { \text{tot} } ( \theta ) & = \ln \mathcal { L } _ { \text{LOS} } + \ln \mathcal { L } _ { \text{PM}, t } + \ln \mathcal { L } _ { \text{PM}, R } + \ln \mathcal { L } _ { \text{VSP1} } \\ & \quad + \ln \mathcal { L } _ { \text{VSP2} } + \ln \mathcal { L } _ { p, \text{LOS} } + \ln \mathcal { L } _ { \text{pos} } + \ln \mathcal { L } _ { \text{pos}, R }, \quad ( 3 9 )$$

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 di ff erent 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,

Table 2. Priors implemented in dynesty for our analysis.

| Parameter    | Units   | Prior Type   | Range                 |
|--------------|---------|--------------|-----------------------|
| e β 0        | none    | flat         | [ - 1 , 1]            |
| e β ∞        | none    | flat         | [ - 1 , 1]            |
| log 10 r t   | pc      | flat         | [0 , 2]               |
| η            | none    | flat         | [0 , 4]               |
| r cen        | pc      | flat         | [0 . 5 , 2 . 5]       |
| log 10 M cen | M ⊙     | flat         | [0 , 6]               |
| log 10 M BH  | M ⊙     | flat         | [0 , 6]               |
| log 10 M ⋆   | M ⊙     | flat         | [5 , 8]               |
| log 10 M p   | M ⊙     | flat         | [0 , 6]               |
| α            | none    | flat         | [ - 8 , - 2]          |
| r 0          | pc      | flat         | [0 , 10]              |
| log 10 B i   | G       | normal       | µ = 8 . 47 σ = 0 . 33 |
| l            | 4.54 pc | flat         | [ - 10 , 10]          |

Notes. To ensure a broad range of values is considered, the variable log 10 Bi [G] is sampled over the flat range [6, 10], inside which the normal prior is applied. See discussion in text for more details on the priors adopted.

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 di ff ers from the results obtained from a separate treatment.

The prior on l is expressed in terms of the core-radius length scale rc ≡ 4 54 pc of Baumgardt &amp; Hilker (2018), which is also . used by Dai et al. (2023). Following Read &amp; Steger (2017), we also work with the symmetrized anisotropy:

$$\widetilde { \beta } & \equiv \frac { \rho } { 2 - \beta } & ( 4 0 ) \overset { \sim \ 1 U \ 1 V } { \text{matically} }$$

for -1 ≤ e β ≤ 1 to e ffi ciently 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). 9 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 su ffi ciently narrow to allow for satisfactory numerical e ffi ciency. Priors in log-space are designed to e ffi ciently 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

9 We found that extending these limits within the physically allowable range that allows for convergence and positivity of the distribution does not greatly a ff ect its overall shape due to the degeneracy in the parameters r 0 and α.

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 e ffi cient 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. 10

## 6. Kinematic fits and constraints on a central mass

Fig. 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 χ ν 2 = 1 61 for our 4-parameter anisotropy profile and a 3-. parameter mass model (stellar mass and Plummer mass and scale radius). 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. 11 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 3 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 4 M and ⊙ ∼ 10 5 M total mass values required for kine-⊙ matically 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 Sec. 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 5 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 3 M ⊙

10 For further information, see:

https://dynesty.readthedocs.io/en/latest/index.html .

11 More precisely, the introduction of MSP and IMBH components leads to a higher χν 2 value, due to the introduction of 4 additional parameters and a negligible e ff ect on the total chi-square statistic, showing no statistically significant preference for these two components based on the additional degrees of freedom that they introduce.

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 Sec. 3. All data, including the photometric profile, have been self-consistently binned with the center from Anderson &amp; Van der Marel (2010). Upper left: Tangential proper motion velocity dispersion. Upper right: Radial proper motion velocity dispersion. Middle left: LOS velocity dispersion. Middle right: 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 Sec. 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 left: 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 Sec. 3). Lower right: The same as the lower left figure, but for the virial shape parameter 2.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-120e6327-6646-48b1-a158-4b13e829cf76.jpg)

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 4 M ⊙ , strongly limiting the kinematic contribution of the more extended MSP component. 12

12 This result can be more dependent on the modeling of the other extended components in the distribution, especially if these exhibit deArticle number, page 9 of 20

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-f0b81172-42e2-470d-9b44-60c0d76cfb25.jpg)

/circledot

/circledot

Fig. 2. 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 | e β | &lt; 0 1 . , 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.

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. 13 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.

generacies with each other. However, this is not the case in our analysis, as a kinematically relevant contribution is disfavored for this component. Further, we found that the precise values the 3 σ bounds in this study can show moderate stochastic variability when comparing multiple independent runs of the fits, while this has only a moderate e ff ect (agreeing at the specified number of significant figures in most cases), we chose the cited values conservatively so that they exclude at least at the 3 σ level, across the multiple (more than 10) runs considered.

13 Note also that while a single-component Plummer sphere or a similar extended mass model can in principle accommodate for a point mass in the limit of vanishing scale radius, it su ff ers the limitation of being a mutually exclusive model in that it does not explore the coexistence between a point mass and an extended one.

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 su ffi ciently elevated within a few pc from the center to favor a significant extended central mass component of ∼ 2 -3 × 10 5 M that is needed in addition to the photometric ⊙ one. Large ( ≳ 10 4 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 su ffi cient 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,

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 × 10 4 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-9c84315d-8436-45f0-9343-38437f4a94c0.jpg)

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 6 M ⊙ , which for a total luminosity of ∼ 1 25 . × 10 6 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 χ ν 2 = 1 13 for the . 4-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.

Fig. 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 una ff ected (cf. Appendix B, where a fit is performed without including MSP accelerations).

The dominant contribution to this e ff ect appears to be attributable to the most distant MSP acceleration data point, which e ff ectively 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 e ff ected by the outermost data point is robust with respect to intrinsic spin-down e ff ects, 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 Sec. 8.2), and by our statistical treatment of intrinsic spin-down modeling which allows quantifying the significance of putative tensions.

For illustrative purposes, we also show in Fig. 3 the e ff ect of concentrating ∼ 15% of the central mass component in the form a 4 × 10 4 M IMBH, which is representative of some of the val-⊙ ues 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. 14 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).

14 Wedetermine this region as that for which our analysis allows a measurably large extremal acceleration that could no be accounted by the extended components alone. This corresponds to the blue, thin inner cusp seen in Fig. 3.

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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-07446c2d-0caf-41fb-8b4c-1c97d6a7d7b8.jpg)

## 7. Analysis of the pulsar distribution

Fig. 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 &amp; 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 &amp; Freire (2014)). These form dynamically in GCs at a rate set, to leading order, by the

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, non-linear, 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 &amp; 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).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-261fda4b-c836-45a2-91e5-15b87f5a8908.jpg)

rate of stellar encounters. 15 Such models predict that the abundance of MSPs in GCs should scale with the (total) stellar encounter rate for the cluster Γ ∝ R dr r 2 ρ⋆ ( r ) 2 /σ ( r ) , where ρ⋆ ( r ) is the 3D density profile of stars, σ ( r ) is the 3D velocity dispersion profile of the stars, with ρ⋆ ( r ) 2 /σ ( r ) being the encounter rate density, which encapsulates the rate of two-body stellar interactions (hence the ρ ⋆ 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 di ff er from the single-star one (e.g. Verbunt &amp; 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 &amp; Freire (2014)), and the e ff ects 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 &amp; 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 intracluster level. In particular, we explore whether the MSP popu-

15 This may also apply more generally to comparable models involving formation via two-body encounters that are not limited to low mass Xray binaries.

lation density scales with the stellar encounter rate density so that ρ p ( r ) ∝ ρ⋆ ( r ) 2 /σ ( 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:

$$\Gamma ( < R ) \in \int _ { - \infty } ^ { \infty } d l \int _ { 0 } ^ { R } d R ^ { \prime } R ^ { \prime } \rho _ { * } ( r ) ^ { 2 } / \sigma ( r ), & & ( 4 1 ) & \text{degenera} \\ & \text{general}$$

where r = √ R ′ 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. Γ ( &lt; ∞ ∼ ) 90 4). For obtain-. ing 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{Number of MSPs} \, ( < R ) = K \, \Gamma ( < R ), & & ( 4 2 ) \quad _ { \substack { \iota \iota \iota \iota \iota \iota \\ \ln \alpha \ r a c e a } } ^ { ( \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \iota \jmath }$$

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 &amp; 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 4 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 di ff erent behavior, with a distinct non-linear dependence followed in both cases. This further illustrates the di ff erences between these distributions, with the linear scaling behavior being a unique characteristic of the MSP population.

To the knowledge of the authors, while such scaling relations have been explored at the level of total encounter rates and MSPcounts for GC populations, this constitutes the first analysis where a similar relationship is validated as a function of radius within a single GC, thus motivating a functional form for the profile of the pulsar distribution within the cluster.

It is also interesting to note that, following the discussion of Sec. 8.2, since significant mass segregation for lighter stellar remnants (which MSPs are part of) is not expected for ω Cen , a di ff erent 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. 16

## 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 di ff erences 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 e ff ects with anisotropic mass components (e.g. Zocchi et al. (2017); Read &amp; Steger (2017)). As noted in Sec. 8.2, degenerate e ff ects 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 selfconsistent 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 &amp; 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 di ff erent datasets across multiple studies, including those of Anderson &amp; 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 &amp; 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 di ff erent 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 &amp; 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 e ff ects 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)

16 Cf. Chen et al. (2023), where the idealness of ω Cen for probing the MSP distribution based on its lack of mass segregation was also independently argued as a promising direction for future studies.

center. This shows that the essential conclusions of our analysis are robust with respect to the choice of centers considered, with the main e ff ect 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 di ff erent distance dependencies of PM and LOS velocity components and adds further support to the distance determinations in Watkins et al. (2015b), which exploit this e ff ect. This value also shows agreement with the distance determinations of the Harris (1996) catalog (2010 edition) and the more recent ones from Baumgardt &amp; Vasiliev (2021).

Van der Marel &amp; Anderson (2010) reported central values for the central and outer anisotropy parameters of β 0 = 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 di ff erences 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 &amp; Anderson (2010), would be the photometric profile determinations. These authors found that di ff erent assumptions about the morphology of the photometric component led to di ff erences 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. 17

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 e ff ects 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). 18 After making these

17 One such example we considered was a generic three-component Plummer model, which we found to produce significantly inferior fits to the photometry, despite having a larger number of free parameters.

18 Note that the definition of anisotropy is not necessarily equivalent between axisymmetric and spherical models (cf. Eq. 6 with Equation 3 of Evans et al. (2022)), nonetheless, for a system that is close to isotropic and spherical this should be approximately the case, while also allowing

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 e ff ect 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 di ff erences due to the photometric profiles adopted which may not be intrinsic to these e ff ects. 19 This is instructive to compare with D'Souza &amp; Rix (2013), who concluded that flattened models with rotation led to mass estimates for the cluster that were 10% lower, with a 4 2% di . ff erence in the estimate being attributable to rotation. We also note the discussion of Van der Marel &amp; Anderson (2010), where it is argued that the e ff ects 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 e ff ect for the purposes of this study.

We have also considered the e ff ects of di ff erent mass parametrizations for our dark component, finding that the essential conclusions of our analysis remain unchanged. This is explored in Appendix D.

Lastly, 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 e ff ects 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 e ff ect on the central values, also favoring a central dark component over an IMBH. The main di ff erence 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 Sec. 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 present-day GCs (Kremer

for di ff erences in these definitions to be reabsorbed as di ff erent values for the anisotropy profiles.

19 We did, however, perform an earlier analysis based on the independent dataset of Trager et al. (1995) for the photometry, finding consistency with the main results presented in this work.

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 &amp; 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 stellarmass black holes, with total mass 1 82 . + 0 05 . -0 06 . × 10 5 M ⊙ , 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, 20 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 8 200 M , ⊙ based on the observation of centrally located fast-moving stars, in apparent tension with our 3 σ bound of 6 × 10 3 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 Sec. 8.2).

20 We note that, while the modeling of Dickson et al. (2024) indicates a more concentrated distribution for the remnants below this radius which is not reproduced by the Plummer model, our inclusion of a point mass has the ability to compensate for this should a higher enclosed mass be favored for this region. The fact that it does not do so at a statistically significant level indicates that such feature, while physically expected from the remnant models of Dickson et al. (2024), is likely not relevant for the kinematics considered in our analysis.

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 4 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 3 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 4 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 &amp; 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, 21 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 contents 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 which are not present in previous works: One is that a two-component model with a significantly extended central mass distribution (with scale radius 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 other result, which is a consequence of the previous one, 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

21 Note that this does not a priori imply the absence of an IMBH, as our derived bounds are still compatible with lighter IMBHs.

a 6 × 10 3 M 3 ⊙ σ CL upper limit on its mass. This altogether favors 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 4 M ⊙ , was found not to contribute meaningfully to the dynamical mass distribution of ω Cen , we found that the MSP timing constraints do play a significant role in the kinematics, favoring a central distribution that is ∼ 20% more massive and extended.
- -We found that the MSP density profile goes as:

$$\rho _ { p } ( r ) \in \rho _ { * } ( r ) ^ { \gamma } / \sigma ( r ),$$

with γ = 1 9 . ± 0 3 . , consistent with models in which MSPs originate from stellar encounters. It is a natural expectation from this scenario that γ ∼ 2 . We found, further, that the stars and X-ray sources in ω Cen are more radially extended than the MSPs and follow a similar density profile to one another. Our analysis leads to a significant novelty: that the encounter rate-MSP abundance correlation can be extended at 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 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 shows the potential and promise of combining velocity dispersion data with pulsar accelerations to probe the inner mass distribution of GCs. Future observations may allow probing whether extremal accelerations occur for central MSPs, a 'smoking-gun' signature of an IMBH and a disfavored scenario for ω Cen in our analysis.
- -With the number of MSP discoveries in GCs having recently doubled over a 5-year period amid the advent of radio surveys with unprecedented sensitivity, such as MeerKAT, FAST, and notably the upcoming SKA, it is foreseeable that the methodology we introduce will be increasingly relevant for the understanding of the kinematics of these systems.

## References

Abbate, F., Possenti, A., Colpi, M., &amp; Spera, M. 2019, ApJ, 884, L9 Abbate, F., Possenti, A., Ridolfi, A., et al. 2018, Mon. Not. Roy. Astron. Soc., 481, 627 Anderson, J. &amp; Van der Marel, R. P. 2010, The Astrophysical Journal, 710, 1032 Anderson, S. B. 1993, A study of recycled pulsars in globular clusters (California Institute of Technology) Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2023, MNRAS, 526, 429 Baes, M. &amp; van Hese, E. 2007, A&amp;A, 471, 419 Bagchi, M., Lorimer, D. R., &amp; Chennamangalam, J. 2011, MNRAS, 418, 477 Bahramian, A., Heinke, C. O., Sivako ff , G. R., &amp; Gladstone, J. C. 2013, ApJ, 766, 136 Battaglia, G., Helmi, A., &amp; Breddels, M. 2013, New Astronomy Reviews, 57, 52 Baumgardt, H. 2017, MNRAS, 464, 2174 Baumgardt, H., He, C., Sweet, S. M., et al. 2019, Monthly Notices of the Royal Astronomical Society, 488, 5340 Baumgardt, H. &amp; Hilker, M. 2018, Monthly Notices of the Royal Astronomical Society, 478, 1520 Baumgardt, H. &amp; Vasiliev, E. 2021, MNRAS, 505, 5957 Bekki, K. &amp; Freeman, K. C. 2003, Mon. Not. Roy. Astron. Soc., 346, L11

- Bellini, A., Anderson, J., Bedin, L. R., et al. 2017, ApJ, 842, 6
- Bellini, A., Anderson, J., van der Marel, R. P., et al. 2014, ApJ, 797, 115 Binney, J. &amp; Mamon, G. A. 1982, MNRAS, 200, 361 Binney, J. &amp; Tremaine, S. 2008, Galactic Dynamics: Second Edition Brown, A. M., Massey, R., Lacroix, T., et al. 2019 [ arXiv:1907.08564 ] Cappellari, M. 2008, MNRAS, 390, 71 Chen, W. et al. 2023, Mon. Not. Roy. Astron. Soc., 520, 3847 Chon, S. &amp; Omukai, K. 2020, MNRAS, 494, 2851 Collins, M. L. M., Read, J. I., Ibata, R. A., et al. 2021, MNRAS, 505, 5686 Dai, S., Johnston, S., Kerr, M., et al. 2023, Mon. Not. Roy. Astron. Soc., 521, 2616

Dai, S., Johnston, S., Kerr, M., et al. 2020, Astrophys. J. Lett., 888, L18 de Boer, T. J. L., Gieles, M., Balbinot, E., et al. 2019, MNRAS, 485, 4906 de Lorenzi, F., Gerhard, O., Coccato, L., et al. 2009, MNRAS, 395, 76 Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2014, MNRAS, 441, 2986 Dickson, N., Hénault-Brunet, V., Baumgardt, H., Gieles, M., &amp; Smith, P. J. 2023, MNRAS, 522, 5320

Dickson, N., Smith, P. J., Hénault-Brunet, V., Gieles, M., &amp; Baumgardt, H. 2024, MNRAS, 529, 331

Dinescu, D. I., Girard, T. M., &amp; van Altena, W. F. 1999, AJ, 117, 1792 D'Souza, R. &amp; Rix, H.-W. 2013, MNRAS, 429, 1887

- Evans, A. J., Strigari, L. E., &amp; Zivick, P. 2022, Mon. Not. Roy. Astron. Soc., 511, 4251

Feroz, F., Hobson, M. P., &amp; Bridges, M. 2009, MNRAS, 398, 1601

Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24

Freire, P. C. C. et al. 2017, Mon. Not. Roy. Astron. Soc., 471, 857

Freundlich, J., Jiang, F., Dekel, A., et al. 2020, MNRAS, 499, 2912

Fritz, T. K., Chatzopoulos, S., Gerhard, O., et al. 2016, ApJ, 821, 44

- Fujii, M. S., Wang, L., Tanikawa, A., Hirai, Y., &amp; Saitoh, T. R. 2024, arXiv eprints, arXiv:2406.06772

Genina, A., Read, J. I., Frenk, C. S., et al. 2020, MNRAS, 498, 144 Gieles, M. &amp; Zocchi, A. 2015, MNRAS, 454, 576

Gray, E. I., Read, J. I., Taylor, E., et al. 2024, arXiv e-prints, arXiv:2405.19286 Grindlay, J. E., Camilo, F., Heinke, C. O., et al. 2002, Astrophys. J., 581, 470 Handley, W. J., Hobson, M. P., &amp; Lasenby, A. N. 2015a, MNRAS, 450, L61 Handley, W. J., Hobson, M. P., &amp; Lasenby, A. N. 2015b, MNRAS, 453, 4384 Harris, C. R., Millman, K. J., van der Walt, S. J., et al.

- 2020, Nature, 585, 357-362
- Harris, W. E. 1996, AJ, 112, 1487
- Heinke, C. O., Grindlay, J. E., Edmonds, P. D., et al. 2005, ApJ, 625, 796

Helmi, A., Van Leeuwen, F., McMillan, P., et al. 2018, Astronomy &amp; astro- physics, 616, A12

Hénault-Brunet, V., Gieles, M., Strader, J., et al. 2020, MNRAS, 491, 113 Henleywillis, S., Cool, A. M., Haggard, D., et al. 2018, MNRAS, 479, 2834 Hernquist, L. 1990, ApJ, 356, 359

- Higson, E., Handley, W., Hobson, M., &amp; Lasenby, A. 2019, Statistics and Computing, 29, 891
- Hunter, J. D. 2007, Computing in Science &amp; Engineering, 9, 90
- Häberle, M., Neumayer, N., Seth, A., et al. 2024, Fast-moving stars around an intermediate-mass black hole in Omega Centauri

Ibata, R. A., Bellazzini, M., Malhan, K., Martin, N., &amp; Bianchini, P. 2019, Nature

Astronomy, 3, 667

- Jalali, B., Baumgardt, H., Kissler-Patig, M., et al. 2012, Astronomy &amp; Astrophysics, 538, A19
- Jeans, J. H. 1922, MNRAS, 82, 122
- Johnson, C. I., Dupree, A. K., Mateo, M., et al. 2020, The Astronomical Journal, 159, 254
- Júlio, M. P., Brinchmann, J., Zoutendijk, S. L., et al. 2023 [ arXiv:2307.13035 ]
- Kızıltan, B., Baumgardt, H., &amp; Loeb, A. 2017a, Nature, 542, 203
- Kızıltan, B., Baumgardt, H., &amp; Loeb, A. 2017b, Nature, 545, 510
- Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in IOS Press, 87-90

Koposov, S., Speagle, J., Barbary, K., et al. 2022, joshspeagle dynesty: v1.2.3

/

Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., &amp; Rasio, F. A. 2020, in

Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia,

M. Davies, A. Sills, &amp; E. Vesperini, Vol. 351, 357-366

Łokas, E. L. &amp; Mamon, G. A. 2003, MNRAS, 343, 401

Lugger, P. M., Cohn, H. N., &amp; Grindlay, J. E. 1995, ApJ, 439, 191

Mamon, G. A. &amp; Łokas, E. L. 2005, MNRAS, 363, 705

Manchester, R. N., Hobbs, G. B., Teoh, A., &amp; Hobbs, M. 2005, VizieR Online Data Catalog, VII 245 /

Massari, D., Koppelman, H. H., &amp; Helmi, A. 2019, Astronomy &amp; Astrophysics, 630, L4

- Merrifield, M. R. &amp; Kent, S. M. 1990, AJ, 99, 1548

Merritt, D. 1985, MNRAS, 214, 25P

- Myeong, G., Vasiliev, E., Iorio, G., Evans, N., &amp; Belokurov, V. 2019, Monthly Notices of the Royal Astronomical Society, 488, 1235
- Neal, R. M. 2003, The Annals of Statistics, 31, 705
- Nice, D. J. &amp; Taylor, J. H. 1995, ApJ, 441, 429
- Noyola, E., Gebhardt, K., &amp; Bergmann, M. 2008, Astrophys. J., 676, 1008
- Noyola, E., Gebhardt, K., Kissler-Patig, M., et al. 2010, ApJ, 719, L60

Osipkov, L. P. 1979, Pisma v Astronomicheskii Zhurnal, 5, 77

- Padmanabh, P. V., Ransom, S. M., Freire, P. C. C., et al. 2024, arXiv e-prints, arXiv:2403.17799

Pechetti, R., Kamann, S., Krajnovi´ c, D., et al. 2024, MNRAS, 528, 4941

- Perera, B. B. P., Stappers, B. W., Lyne, A. G., et al. 2017a, MNRAS, 471, 1258
- Perera, B. B. P., Stappers, B. W., Lyne, A. G., et al. 2017b, Mon. Not. Roy. Astron. Soc., 468, 2114
- Phinney, E. S. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 50, Structure and Dynamics of Globular Clusters, ed. S. G. Djorgovski &amp;G. Meylan, 141

Plummer, H. C. 1911, MNRAS, 71, 460

- Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., &amp; McMillan, S. L. W. 2004, Nature, 428, 724
- Prager, B. J., Ransom, S. M., Freire, P. C. C., et al. 2017, Astrophys. J., 845, 148

Read, J. I. &amp; Steger, P. 2017, MNRAS, 471, 4541

Read, J. I., Walker, M. G., &amp; Steger, P. 2018, MNRAS, 481, 860

- Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130
- Reijns, R. A., Seitzer, P., Arnold, R., et al. 2006, A&amp;A, 445, 503
- Reynoso-Cordova, J., Burgueño, O., Geringer-Sameth, A., et al. 2021, JCAP, 02, 010

Richardson, T. &amp; Fairbairn, M. 2014, MNRAS, 441, 1584

Ridolfi, A., Gautam, T., Freire, P. C. C., et al. 2021, MNRAS, 504, 1407

Shklovskii, I. 1970, Soviet Astronomy, Vol. 13, p. 562, 13, 562

Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, &amp; U. V. Toussaint, 395-405

Skilling, J. 2006, Bayesian Analysis, 1, 833

- Smith, P. J., Hénault-Brunet, V., Dickson, N., Gieles, M., &amp; Baumgardt, H. 2024, arXiv e-prints, arXiv:2407.06274
- Sollima, A., Bellazzini, M., Smart, R. L., et al. 2009, MNRAS, 396, 2183
- Speagle, J. S. 2020, MNRAS, 493, 3132
- Strigari, L. E., Bullock, J. S., &amp; Kaplinghat, M. 2007, ApJ, 657, L1
- Trager, S. C., King, I. R., &amp; Djorgovski, S. 1995, AJ, 109, 218
- van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., &amp; de Zeeuw, P. T. 2006, A&amp;A, 445, 513
- van der Marel, R. P. 1994, MNRAS, 270, 271

Van der Marel, R. P. &amp; Anderson, J. 2010, The Astrophysical Journal, 710, 1063

- Vasiliev, E. &amp; Baumgardt, H. 2021, MNRAS, 505, 5978

Verbunt, F. &amp; Freire, P. C. C. 2014, A&amp;A, 561, A11

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

- Vleeschower, L., Corongiu, A., Stappers, B. W., et al. 2024, MNRAS[ arXiv:2403.12137 ]

Wang, L.-C. &amp; Xie, Y. 2021, Research in Astronomy and Astrophysics, 21, 270

- Watkins, L. L., van de Ven, G., den Brok, M., &amp; van den Bosch, R. C. E. 2013, MNRAS, 436, 2598
- Watkins, L. L., van der Marel, R. P., Bellini, A., &amp; Anderson, J. 2015a, ApJ, 803, 29
- Watkins, L. L., van der Marel, R. P., Bellini, A., &amp; Anderson, J. 2015b, ApJ, 812, 149
- Wilkinson, M. I., Kleyna, J., Evans, N. W., &amp; Gilmore, G. 2002, Mon. Not. Roy. Astron. Soc., 330, 778
- Wirth, H., Bekki, K., &amp; Hayashi, K. 2020, Monthly Notices of the Royal Astronomical Society: Letters, 496, L70

Xie, Y. &amp; Wang, L.-C. 2020, Research in Astronomy and Astrophysics, 20, 191 Xie, Y., Yin, D., Wang, L., et al. 2024, MNRAS, 527, 7743

Yin, D., Zhang, L.-y., Qian, L., et al. 2024, arXiv e-prints, arXiv:2405.18228

Zentner, A., Dandavate, S., Slone, O., &amp; Lisanti, M. 2022, JCAP, 07, 031

Zhang, L., Freire, P. C. C., Ridolfi, A., et al. 2023, ApJS, 269, 56

Zhao, H. 1996, MNRAS, 278, 488

Zhao, J. &amp; Heinke, C. O. 2022, MNRAS, 511, 5964

Zhao, J. &amp; Heinke, C. O. 2023, MNRAS, 526, 2736

Zocchi, A., Gieles, M., &amp; Hénault-Brunet, V. 2017, MNRAS, 468, 4429

Zocchi, A., Gieles, M., &amp; Hénault-Brunet, V. 2019, MNRAS, 482, 4713

## 10. 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énaultBrunet, 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 highperformance 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 DrlicaWagner 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. AB, FC and JMC acknowledge funding received from the European Union through the grant 'UNDARK' of the Widening participation and spreading excellence programme (project number 101159929). AB and JMC acknowledge support from the MICINN through the grant 'DarkMaps' PID2022-142142NB-I00. JIR would like to acknowledge support from STFC grants ST Y002865 1 and / / ST Y002857 1. This work has made use of the following soft-/ / ware packages: dynesty (Speagle 2020; Koposov et al. 2022), GravSphere (Read &amp; Steger 2017; Read et al. 2018; Genina et al. 2020; Collins et al. 2021), corner.py (ForemanMackey 2016), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Jupyter Notebook (Kluyver et al. 2016).

## Appendix A: Posterior distributions of the anisotropy parameters

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-8ea4e625-903e-4618-8dba-5ff2197eab2d.jpg)

∞

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). 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

## Appendix C: Use of different centers

## Appendix D: Use of different dark components

Figs. D.1 and D.2 explore the influence of di ff erences 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. Note that the scale radius (which is mathematically identical to r 0 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

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-7d6fec58-db0a-45d4-8b7a-697415aefbfe.jpg)

/circledot

/circledot

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 di ff erences.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-7ae4dca2-eb05-4dd6-9447-43505bc17abf.jpg)

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 di ff erent posterior distribution when accelerations are excluded. This, however, has a negligible e ff ect on the remainder of the mass components in the absence of accelerations to constrain their kinematics.

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 di ff erences. 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-c41d9065-6109-43f3-94b3-1c969e5fc7bb.jpg)

† Marginal di ff erences 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 di ff erences (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.

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).

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-c5809cd3-e9b5-40e0-979f-f6979cb9e54b.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-c630f896-b49d-4de6-a711-64dc6fd88210.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-845aa8d5-911a-45cc-89e6-3231b9837968.jpg)

Fig. C.2. Combined velocity dispersion profiles for the three components using the Anderson &amp; 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.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-1422b70a-ae7b-475c-a81f-a63a0b9a50f1.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00939v5-figure-363326c7-17e6-4811-945e-c6d63d5b9df9.jpg)

/circledot

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}, g } ( r ) = \frac { M _ { \text{cen} } } { ( 2 \pi r _ { g } ^ { 2 } ) ^ { 3 / 2 } } \exp \left ( - \frac { r ^ { 2 } } { 2 r _ { g } ^ { 2 } } \right ).$$

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 5 M and the absence of a kine-⊙ matically 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.
Andrés Bañares-Hernández Francesca Calore Jorge Martin Camalich Justin I. Read
New constraints on the central mass contents of Omega Centauri from combined stellar kinematics and pulsar timing
We perform a combined analysis of stellar kinematics and line-of-sight
accelerations of millisecond pulsars to investigate the mass contents of Omega
Centauri. We consider multiple mass components: the visible photometric
distribution, a central cluster of dark remnants, an intermediate-mass black
hole, and a distribution of intermediate extension traced by the pulsars. By
self-consistently incorporating multiple independent datasets, including the
effects of different centers, we obtain significant constraints on these mass
distributions. Our results strongly favor an extended central mass of ~ 2 - 3 x
10^5 Msun, emulating a cluster of heavy stellar remnants, over an IMBH, with a
3-sigma upper limit of 6 x 10^3 Msun on its mass. Pulsar timing observations
provide significant constraints, favoring a central mass distribution that is ~
20% more massive and extended. Notably, we observe a clear linear scaling
between the pulsar distribution and stellar encounter rates, showing excellent
agreement with the expectation from millisecond pulsar formation models and
motivating a profile for their distribution, providing the first validation of
its kind where positions are linked to their place of formation in globular
clusters. Our findings highlight the effectiveness of this novel methodology in
exploring the structure of star clusters, setting a promising precedent amid
the advent of the rapidly growing number of observations and discoveries
currently being experienced in this field.