We perform a combined analysis of stellar kinematics and lineofsight 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 intermediatemass black hole and a distribution of intermediate extension traced by the pulsars. By selfconsistently 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 105 Msun emulating a cluster of heavy stellar remnants over an IMBH with a 3sigma upper limit of 6 x 103 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.
['Astrophysics - Astrophysics of Galaxies']
New constraints on the central mass contents of Omega Centauri from combined stellar kinematics and pulsar timing
{'No Header': '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  \n- 1 Instituto de Astrofísica de Canarias, La Laguna, Tenerife, E-38200, Spain\n- 2 Departamento de Astrofísica, Universidad de La Laguna\n- 3 CNRS, Laboratoire d'Annecy-le-Vieux de Physique Théorique, 74940, Annecy, France\n- 4 CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland\n- 5 Department of Physics, University of Surrey, Guildford, GU2 7XH, United Kingdom  \nNovember 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.  \nMethods. 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 .  \nResults. 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.  \nConclusions. 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.  \nKey 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 & 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  \nbeen seen in several massive Milky Way GCs, including ω Cen (Gray et al. 2024).  \nRecent 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).  \nAlongside the possible presence of dark matter, ω Cen is interesting as a possible location for the formation of black hole seeds via runaway stellar collisions (e.g. Portegies Zwart et al. 2004) or the formation and death of a supermassive star (e.g. Chon & Omukai 2020). Several authors (Noyola et al. 2008; Noyola et al. 2010; Jalali et al. 2012; Baumgardt 2017) found evidence of a ∼ 3 -5 × 10 4 M ⊙ intermediate-mass black hole  \n(IMBH) at the center of ω Cen from the consideration of stellar kinematic data (using either Jeans models or N-body simulations). However, Van der Marel & Anderson (2010), via a 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.  \nIn 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 & Wang 2020; Ridolfi et al. 2021; Wang & Xie 2021; Zhang et al. 2023; Dai et al. 2023; Padmanabh et al. 2024; Xie et al. 2024; Vleeschower et al. 2024). We also refer the interested reader to earlier studies where this methodology was explored, such as in Phinney (1993) and Anderson (1993).  \nIn 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 particular, 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.  \nWith 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 & 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  \nAn 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)).  \nThis 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 & Steger (2017); Genina et al. (2020); Collins et al. (2021) and Júlio et al. (2023).  \nIn 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.  \nWe 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:  \nd f d t = ∂ f ∂ t + ∇ x f · v - ∇ v f · ∇ x Φ = 0 , (1)  \nwhere Φ is the gravitational potential resulting from all mass components in the system, and satisfying the Poisson equation:  \n∇ 2 x Φ = 4 π G ρ, (2)  \nwhere ρ is the total mass density of the system.  \nAlthough 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 & 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 j , one obtains a set of three wellknown equations: The Jeans equations (Jeans 1922), which can be expressed as:  \n3 X i = 1 ∂ν⋆ ⟨ v i v j ⟩ ∂ x i + ν⋆ ∂ Φ ∂ x j = 0 , (3)  \nimplement this within the distribution function modeling code limepy (Gieles & Zocchi 2015). They then apply their method to the 47 Tuc and Terzan 5 GCs, whereas we focus here on ω Cen.  \nwhere ν⋆ ( x ) = R d 3 v f is the stellar (tracer particle) number density and ⟨ v i v j ⟩ = R d 3 vv i v j f correspond to the two-point velocity moments, with the brackets \' ⟨⟩ \' generally denoting an integral with the distribution function over velocity space.  \nFor 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):  \nσ 2 r / t = ⟨ v 2 r / t ⟩ - ⟨ v r / t ⟩ 2 (4)  \nto give the spherical Jeans equation:  \n1 ν⋆ ∂ν⋆σ 2 r ∂ r + 2 βσ 2 r r = -GM r 2 , (5)  \nwhere  \nβ ≡ 1 -σ 2 t σ 2 r (6)  \nis 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 .  \nEq. (5) can be solved for the quantity ν⋆σ 2 r yielding the radial velocity dispersion (Mamon & Łokas 2005; van der Marel 1994):  \nσ 2 r ( r ) = 1 ν⋆ ( r ) g ( r ) Z ∞ r GM ( r \' ) ν⋆ ( r \' ) r \' 2 g ( r \' ) dr \' , (7)  \nwhere  \ng ( r ) ≡ exp 2 Z β ( r ) r dr ! . (8)  \nThe LOS projection of σ 2 r , which is a commonly used observable (e.g. Strigari et al. (2007)), is then given by (Binney & Mamon 1982)  \nσ 2 LOS ( R ) = 2 Σ ⋆ ( R ) Z ∞ R GLYPH<18> 1 -R 2 r 2 β ( r ) GLYPH<19> ν⋆ ( r ) σ 2 r ( r ) r √ r 2 -R 2 dr , (9)  \nwhere Σ ∗ ( 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.  \nHaving 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 & Kent 1990; Wilkinson et al. 2002; Łokas & Mamon 2003; de Lorenzi et al. 2009; Read & Steger 2017). This problem is a manifestation of the fact that the Jeans equations alone are not a closed system, meaning that they cannot unambiguously constrain these two quantities without more information being provided.  \nFollowing, for example, Watkins et al. (2013); Van der Marel & Anderson (2010), in order to resolve the mass-anisotropy degeneracy, one can obtain expressions for the projected proper motions along the projected radius R and the tangential component (respectively), these are given by:  \nσ 2 PM , R ( R ) = 2 Σ ⋆ ( R ) Z ∞ R GLYPH<18> 1 -β ( r ) + R 2 r 2 β ( r ) GLYPH<19> ν⋆ ( r ) σ 2 r ( r ) r √ r 2 -R 2 dr (10)  \nand  \nσ 2 PM , t ( R ) = 2 Σ ⋆ ( R ) Z ∞ R GLYPH<18> 1 -β ( r ) GLYPH<19> ν⋆ ( r ) σ 2 r ( r ) r √ r 2 -R 2 dr . (11)  \nFurthermore, following Merrifield & Kent (1990); Richardson & Fairbairn (2014), to help eliminate the ρ -β degeneracy, one can integrate with higher moments of velocity in the collisionless Boltzmann equation, Eq. (1), to obtain two independent observables known as virial shape parameters (VSPs), which are given by  \nVSP1 = 2 5 Z ∞ 0 GM ν⋆ (5 -2 β ) σ 2 r rdr = Z ∞ 0 Σ ⋆ ⟨ v 4 LOS ⟩ RdR (12)  \nand  \nVSP2 = 4 35 Z ∞ 0 GM ν⋆ (7 -6 β ) σ 2 r r 3 dr = Z ∞ 0 Σ ⋆ ⟨ v 4 LOS ⟩ R 3 dR , (13)  \nwhere one can obtain data from observations corresponding to the RHS of Eqs. (12) and (13) that can therefore constrain β.  \nFor 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:  \nν⋆ ( r ) = ρ c ( r / rc ) γ (1 + ( r / rc ) α ) ( γ -β ) /α , (14)  \nwhere 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.  \nThe projected tracer surface density, is then given by:  \nΣ ⋆ ( R ) = Z ∞ -∞ ν⋆ ( r ( l )) dl , (15)  \nwhere, 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.  \nTo 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.  \nFor 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  \nM cen( r ) = M cen r 3 r 3 cen 1 + r 2 r 2 cen ! -3 / 2 , (16)  \nwhere 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 .  \nLastly, 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  \nMp ( r ) = 4 Mp 3 √ π Γ h (1 -α ) / 2 i Γ h -( α/ 2 + 1) i r 3 r 3 0 2 F 1 " 3 2 , 1 -α 2 ; 5 2 ; -r 2 r 2 0 # , (17)  \nwhere 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 α < -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:  \nM ( r ) = M ⋆ ( r ) + M cen( r ) + M BH + Mp ( r ) . (18)  \nWhile 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  \ncase, a more concentrated distribution of pulsars could be explained by their mechanism of formation rather than mass segregation (see Sec. 7).  \nThere 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  \nWe use the anisotropy profile assumed in GravSphere following the generalized form from Baes & van Hese (2007), which is an extension of the Osipkov-Merritt anisotropy profile (Osipkov 1979; Merritt 1985), and is given by  \nβ ( r ) = β 0 + ( β ∞ -β 0) 1 1 + GLYPH<16> rt r GLYPH<17> η , (19)  \nwhere β 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 & 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 & 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.  \nThe LOS velocity and PM data are then binned using the binulator , as described in Collins et al. (2021). We use 100  \nlogarithmically 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:  \n˙ P P ! obs = a LOS + a S + a g c + ˙ P P ! int , (20)  \nwhere 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.  \nag 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 & Taylor (1995) yielding  \nag = -cos( b ) Θ 2 0 R 0 ! cos( ℓ ) + ϑ sin 2 ( ℓ ) + ϑ 2 ! ms -2 , (21)  \nwhere ϑ ≡ ( 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 chosen for consistency with the values reported in Bellini et al. (2014); Watkins et al. (2015a,b), and was subsequently validated throughout our analysis (see 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  \nFor 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.  \naS 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)  \naS = 3 . 78 × 10 -12 d 5 . 2 kpc ! µ T mas yr -1 ! 2 ms -2 , (22)  \nwhere µ 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).  \na 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  \na LOS( r , l ) = -G l r ! M ( r ) r 2 , (23)  \nwhere 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.  \nWhile 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)  \na int ≡ c ˙ P P ! int ≈ 7 . 96 × 10 -10 B 2 × 10 8 G ! 2 2ms P obs ! 2 ms -2 , (24)  \nwhere 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.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:  \nln L GLYPH<0> θ GLYPH<1> = -1 2 X y χ 2 y , (25)  \nanalysis 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.  \nwhere we have the general chi-square variable for an observable y defined as  \nχ 2 y ≡ X i GLYPH<20> yi , obs -yi GLYPH<0> θ GLYPH<1> GLYPH<21> 2 δ y 2 i , (26)  \nand observables are taken to be normally distributed with median values { yi , obs } and corresponding standard errors { δ yi } , and { yi GLYPH<0> θ GLYPH<1> } 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.  \nTo allow for an e ff ective implementation and with the above discussion in mind, we include in the fit the quantity:  \na obs ≡ c ˙ P P ! obs -ag -aS = a LOS + c ˙ P P ! int , (27)  \nwhere 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.  \nThe term to be added to the log-likelihood function takes the form:  \nln L p , LOS( θ M ; { li } , { Bi } ) = -1 2 Np X i ( a obs , i -a ( θ M ; li , Bi )) 2 δ a 2 i , (28)  \nwhere a ( θ M ; li , 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 { li } 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.  \nTo perform the fit, we marginalize over the parameters { Bi } and { li } . The li 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  \npotential 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 deviation of 0.33 dex, which we adopt as a prior. This allows for an ample range which captures some of the uncertainty in the intrinsic spin-down, while also being representative of the values expected for MSP populations.  \na 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 { li } to determine whether the bound is satisfied.  \nTo 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):  \nn ( R ) = n 0 " 1 + GLYPH<18> R r 0 GLYPH<19> 2 # α/ 2 , (29)  \nwhere 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.)  \nAs seen in Grindlay et al. (2002), the 3D density has the form:  \nn ( r ) = f ( α, r 0 , n 0) " 1 + GLYPH<18> r r 0 GLYPH<19> 2 # ( α -1) / 2 , (30)  \nwhere 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  \nn ( R ) = Z ∞ -∞ n √ R 2 + l 2 ! dl , (31)  \nand one can verify that this expression yields Eq. (29) with  \nf ( α, r 0 , n 0) = Γ [(1 -α ) / 2] √ π Γ ( -α/ 2) n 0 r 0 , (32)  \nTable 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).  \nNotes. Error propagation follows the approach outlined in Sec. 4.  \nwhere Γ denotes the Euler gamma function.  \nTo 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:  \nP ( l , R | α, r 0) = P ( l | R , α, r 0) P ( R | α, r 0) , (33)  \nwhere:  \nP ( l | R , α, r 0) dl = n √ R 2 + l 2 ! dl R ∞ -∞ n √ R 2 + l \' 2 ! dl \' = Γ [(1 -α ) / 2] √ π Γ ( -α/ 2) 1 r 0 h 1 + GLYPH<16> R 2 + l 2 GLYPH<17> / r 2 0 i ( α -1) / 2 h 1 + ( R / r 0) 2 i α/ 2 dl , (34)  \nand: 8  \nP ( R | α, r 0) dR = Rn ( R ) dR R ∞ 0 R \' n ( R \' ) dR \' = -( α + 2) R r 2 0 " 1 + GLYPH<18> R r 0 GLYPH<19> 2 # α/ 2 dR , (35)  \nwhere 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.  \nThe complete positional probability for a given pulsar is thus given by:  \nP ( l , R | α, r 0) = -α + 2 √ π Γ h (1 -α ) / 2 i Γ h -α/ 2 i R r 3 0 GLYPH<20> 1 + ( l 2 + R 2 ) / r 2 0 GLYPH<21> ( α -1) / 2 , (36)  \nyielding the relevant (i.e., non-constant) contribution to the loglikelihood function for Np pulsars:  \nln L pos( { li } , α, r 0) = α -1 2 Np X i ln(1 + ( R 2 i + l 2 i ) / r 2 0 ) + Np ln Γ h (1 -α ) / 2 i Γ h -α/ 2 i + ln[ -( α + 2)] -3 ln r 0 ! . (37)  \nSince 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 N \' p pulsars):  \nln L pos , R( α, r 0) = α 2 N \' p X i ln " 1 + GLYPH<18> Ri r 0 GLYPH<19> 2 # + N \' p (ln[ -( α + 2)] -2 ln r 0) . (38)  \nNote 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  \nln L tot( θ ) = ln L LOS + ln L PM , t + ln L PM , R + ln L VSP1 + ln L VSP2 + ln L p , LOS + ln L pos + ln L pos , R , (39)  \nwhere θ 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,  \nTable 2. Priors implemented in dynesty for our analysis.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.  \nas 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.  \nThe prior on l is expressed in terms of the core-radius length scale rc ≡ 4 . 54 pc of Baumgardt & Hilker (2018), which is also used by Dai et al. (2023). Following Read & Steger (2017), we also work with the symmetrized anisotropy:  \ne β ≡ β 2 -β (40)  \nfor -1 ≤ e β ≤ 1 to e ffi ciently sample the (infinitely) broad range of possible values of β.  \nTable 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.  \nThe 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  \nposterior 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 3parameter 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.  \nOne 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 kinematically signficant IMBH and MSP profile contributions (respectively) exceed the 3 σ limits of our derived posterior distributions, and are thus self-consistently excluded by the fit. We also found that, when fitting only the MSP component to model extended remnants, the MSP distribution was forced to emulate the concentration of our central mass favored by the kinematics, but failed to reproduce the more extended distribution of the MSPs, which is implicitly accounted for in the position likelihood components. This can be seen when comparing the cumulative distribution functions of the various components, which are examined in Sec. 7 (Fig. 4). This result indicates the need for a more concentrated mass component that is distinct from one traced by the pulsars.  \nOur 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 ⊙  \nFig. 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 & 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.  \n<!-- image -->  \nand 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  \nthe kinematic contribution of the more extended MSP component. 12  \nFig. 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 β | < 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.  \n<!-- image -->  \n/circledot  \n/circledot  \nThis 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  \ngeneracies 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.  \nIMBH, in our analysis we are able to consider them simultaneously in a self-consistent fit.  \nOne 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.  \nThe 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,  \nFig. 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.  \n<!-- image -->  \nwhere the three velocity dispersion components are plotted in the same figure, and is also apparent from the (symmetrized) anisotropy profile in Fig. 2.  \nThe 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.  \nFig. 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  \nmass 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).  \nThe 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.  \nIn 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.  \nFor 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 values in the literature (Noyola et al. 2008; Baumgardt 2017). We observe that, while in both models extremal LOS accelerations are allowable, the larger IMBH model would produce significant extremal accelerations over a larger region within the inner ∼ 0 . 5 pc, comparable to the radius of the circle of influence r inf ∼ of the putative IMBH. A detection of such an extremal acceleration within this region would constitute a 'smoking-gun' signature of such an IMBH. Such scenario is disfavored by stellar kinematics in our analysis, except for the very inner region within ≲ 0 . 15 pc, where a smaller IMBH of at most a few thousand solar masses is still allowable. 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).  \nFig. 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 & 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).  \n<!-- image -->  \nFig. 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.  \n<!-- image -->", '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 & 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).  \nWe now consider what the more concentrated MSP profile that we find here can tell us about MSP formation. MSPs have long been thought to descend from low mass X-ray binaries (e.g. see Bahramian et al. (2013); Verbunt & Freire (2014)). These form dynamically in GCs at a rate set, to leading order, by the  \nrate 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 & Freire 2014), binary population statistics, loss of stars from the cluster (impacted by the cluster escape velocity; Yin et al. 2024), neutron star retention fractions (e.g. Verbunt & Freire (2014)), and the 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 & Heinke 2022).  \nPrevious 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-  \nlation 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.  \nIn 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.  \nFig. 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:  \nΓ ( < R ) ∝ Z ∞ -∞ dl Z R 0 dR ' R ' ρ⋆ ( r ) 2 /σ ( r ) , (41)  \nwhere r = √ R ' 2 + l 2 .  \nEncounter rates are usually expressed in arbitrary units, which in our case have adapted to be consistent with Bahramian et al. (2013) for the total value (i.e. Γ ( < ∞ ) ∼ 90 . 4). For obtaining the constant of proportionality K mapping the encounter rate to the number of MSPs, we perform a least-squares fit (Fig. 5: red line) where we model the cumulative number of MSPs as:  \nNumber of MSPs ( < R ) = K Γ ( < R ) , (42)  \nyielding K ∼ 0 . 25 . This closely matches the central value of the total number of MSPs ( ∼ 23) independently predicted by the phenomenological relation from Zhao & Heinke (2022) (Fig. 5: gray, dashed line). However, the 1 σ upper bound of this relation implies more than double the abundance, suggesting that ω Cen could host a significant number of yet-undiscovered MSPs. One can also use the kinematic component of the MSP distribution to derive a strict upper bound on the MSP abundance assuming the totality of the dynamical mass being due to ∼ 1 . 6 M ⊙ MSP binary systems, obtaining ∼ 4 × 10 4 at 3 σ CL (c.f. Fig. 2).  \nThe 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.  \nOn 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.  \nTo 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.  \nIt 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.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 & 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.  \nNoyola et al. (2010) found a dependence of the inferred central mass with the assumed location of the kinematic center in their analysis, with the center used in Anderson & Van der Marel (2010) leading to a lower estimate of the inferred central mass. In a recent study from Pechetti et al. (2024), employing 3 independent center determinations, it was concluded that the favored center led to inferred stellar kinematics determinations that were consistent with the results from di ff erent datasets across multiple studies, including those of Anderson & Van der Marel (2010) using HST data. However, a significant tension was found with the originally assumed center in Noyola et al. (2008). This is also consistent with the findings of Anderson & Van der Marel (2010), where the Noyola et al. (2008) center was ruled out at high confidence and determined to be susceptible to systematic errors not originally considered in the study of Noyola et al. (2008).  \nWe 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 & 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)  \ncenter. 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.  \nAnother 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 & Vasiliev (2021).  \nVan der Marel & 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.  \nAnother factor which may play a role in central mass determinations as seen, for instance, in Van der Marel & Anderson (2010), would be the photometric profile determinations. These authors found that 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  \nTo 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.  \nAs 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  \nadjustments 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 & 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 & 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.  \nWe 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.  \nLastly, 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  \net 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.  \nDickson et al. (2024) performed a recent analysis of ω Cen among other GCs, in which limepy distribution function models (Gieles & Zocchi 2015) were used to obtain realistic distributions for an ample spectrum of species comprising stellar-mass black holes, white dwarfs, neutron stars and main sequence stars. These were subsequently fitted to velocity dispersion and stellar mass function data following the procedure described in Dickson et al. (2023). A key finding from this study is the presence of a substantial central mass component dominated by 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).  \nPechetti 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.  \nIMBH 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.  \nMany works have suggested that ω Cen is the stripped remnant of a nucleated dwarf (Bekki & Freeman 2003; Wirth et al. 2020; Johnson et al. 2020; Gray et al. 2024). If so, it should have a dark matter halo and a correspondingly higher escape velocity. This could explain the fast-moving stars found in Häberle et al. (2024) without the need for such a massive IMBH, 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:  \n- -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.\n- -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  \na 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.  \n- -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.\n- -We found that the MSP density profile goes as:  \nρ p ( r ) ∝ ρ⋆ ( r ) γ /σ ( r ) ,  \nwith γ = 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.  \n- -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).\n- -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.\n- -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., & Spera, M. 2019, ApJ, 884, L9 Abbate, F., Possenti, A., Ridolfi, A., et al. 2018, Mon. Not. Roy. Astron. Soc., 481, 627 Anderson, J. & 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. & van Hese, E. 2007, A&A, 471, 419 Bagchi, M., Lorimer, D. R., & Chennamangalam, J. 2011, MNRAS, 418, 477 Bahramian, A., Heinke, C. O., Sivako ff , G. R., & Gladstone, J. C. 2013, ApJ, 766, 136 Battaglia, G., Helmi, A., & 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. & Hilker, M. 2018, Monthly Notices of the Royal Astronomical Society, 478, 1520 Baumgardt, H. & Vasiliev, E. 2021, MNRAS, 505, 5957 Bekki, K. & Freeman, K. C. 2003, Mon. Not. Roy. Astron. Soc., 346, L11  \n- Bellini, A., Anderson, J., Bedin, L. R., et al. 2017, ApJ, 842, 6\n- Bellini, A., Anderson, J., van der Marel, R. P., et al. 2014, ApJ, 797, 115 Binney, J. & Mamon, G. A. 1982, MNRAS, 200, 361 Binney, J. & 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. & 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  \nDai, 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., & Smith, P. J. 2023, MNRAS, 522, 5320  \nDickson, N., Smith, P. J., Hénault-Brunet, V., Gieles, M., & Baumgardt, H. 2024, MNRAS, 529, 331  \nDinescu, D. I., Girard, T. M., & van Altena, W. F. 1999, AJ, 117, 1792 D'Souza, R. & Rix, H.-W. 2013, MNRAS, 429, 1887  \n- Evans, A. J., Strigari, L. E., & Zivick, P. 2022, Mon. Not. Roy. Astron. Soc., 511, 4251  \nFeroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601  \nForeman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24  \nFreire, P. C. C. et al. 2017, Mon. Not. Roy. Astron. Soc., 471, 857  \nFreundlich, J., Jiang, F., Dekel, A., et al. 2020, MNRAS, 499, 2912  \nFritz, T. K., Chatzopoulos, S., Gerhard, O., et al. 2016, ApJ, 821, 44  \n- Fujii, M. S., Wang, L., Tanikawa, A., Hirai, Y., & Saitoh, T. R. 2024, arXiv eprints, arXiv:2406.06772  \nGenina, A., Read, J. I., Frenk, C. S., et al. 2020, MNRAS, 498, 144 Gieles, M. & Zocchi, A. 2015, MNRAS, 454, 576  \nGray, 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., & Lasenby, A. N. 2015a, MNRAS, 450, L61 Handley, W. J., Hobson, M. P., & Lasenby, A. N. 2015b, MNRAS, 453, 4384  \n- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357-362\n- Harris, W. E. 1996, AJ, 112, 1487\n- Heinke, C. O., Grindlay, J. E., Edmonds, P. D., et al. 2005, ApJ, 625, 796  \nHelmi, A., Van Leeuwen, F., McMillan, P., et al. 2018, Astronomy & astro-  \nphysics, 616, A12  \nHé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  \n- Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2019, Statistics and Computing, 29, 891\n- Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90\n- Häberle, M., Neumayer, N., Seth, A., et al. 2024, Fast-moving stars around an intermediate-mass black hole in Omega Centauri  \nIbata, R. A., Bellazzini, M., Malhan, K., Martin, N., & Bianchini, P. 2019, Nature  \nAstronomy, 3, 667  \n- Jalali, B., Baumgardt, H., Kissler-Patig, M., et al. 2012, Astronomy & Astrophysics, 538, A19\n- Jeans, J. H. 1922, MNRAS, 82, 122\n- Johnson, C. I., Dupree, A. K., Mateo, M., et al. 2020, The Astronomical Journal, 159, 254\n- Júlio, M. P., Brinchmann, J., Zoutendijk, S. L., et al. 2023 [ arXiv:2307.13035 ]\n- Kızıltan, B., Baumgardt, H., & Loeb, A. 2017a, Nature, 542, 203\n- Kızıltan, B., Baumgardt, H., & Loeb, A. 2017b, Nature, 545, 510\n- Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in IOS Press, 87-90  \nKoposov, S., Speagle, J., Barbary, K., et al. 2022, joshspeagle  \n/  \ndynesty: v1.2.3  \nKremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2020, in  \nStar Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia,  \nM. Davies, A. Sills, & E. Vesperini, Vol. 351, 357-366  \nŁokas, E. L. & Mamon, G. A. 2003, MNRAS, 343, 401  \nLugger, P. M., Cohn, H. N., & Grindlay, J. E. 1995, ApJ, 439, 191  \nMamon, G. A. & Łokas, E. L. 2005, MNRAS, 363, 705  \nManchester, R. N., Hobbs, G. B., Teoh, A., & Hobbs, M. 2005, VizieR Online Data Catalog, VII / 245  \nMassari, D., Koppelman, H. H., & Helmi, A. 2019, Astronomy & Astrophysics, 630, L4  \n- Merrifield, M. R. & Kent, S. M. 1990, AJ, 99, 1548  \nMerritt, D. 1985, MNRAS, 214, 25P  \n- Myeong, G., Vasiliev, E., Iorio, G., Evans, N., & Belokurov, V. 2019, Monthly Notices of the Royal Astronomical Society, 488, 1235\n- Neal, R. M. 2003, The Annals of Statistics, 31, 705\n- Nice, D. J. & Taylor, J. H. 1995, ApJ, 441, 429\n- Noyola, E., Gebhardt, K., & Bergmann, M. 2008, Astrophys. J., 676, 1008\n- Noyola, E., Gebhardt, K., Kissler-Patig, M., et al. 2010, ApJ, 719, L60  \nOsipkov, L. P. 1979, Pisma v Astronomicheskii Zhurnal, 5, 77  \n- Padmanabh, P. V., Ransom, S. M., Freire, P. C. C., et al. 2024, arXiv e-prints, arXiv:2403.17799  \nPechetti, R., Kamann, S., Krajnovi'c, D., et al. 2024, MNRAS, 528, 4941  \n- Perera, B. B. P., Stappers, B. W., Lyne, A. G., et al. 2017a, MNRAS, 471, 1258\n- Perera, B. B. P., Stappers, B. W., Lyne, A. G., et al. 2017b, Mon. Not. Roy. Astron. Soc., 468, 2114\n- Phinney, E. S. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 50, Structure and Dynamics of Globular Clusters, ed. S. G. Djorgovski &G. Meylan, 141  \nPlummer, H. C. 1911, MNRAS, 71, 460  \n- Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W. 2004, Nature, 428, 724\n- Prager, B. J., Ransom, S. M., Freire, P. C. C., et al. 2017, Astrophys. J., 845, 148  \nRead, J. I. & Steger, P. 2017, MNRAS, 471, 4541  \nRead, J. I., Walker, M. G., & Steger, P. 2018, MNRAS, 481, 860  \n- Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130\n- Reijns, R. A., Seitzer, P., Arnold, R., et al. 2006, A&A, 445, 503\n- Reynoso-Cordova, J., Burgueño, O., Geringer-Sameth, A., et al. 2021, JCAP, 02, 010  \nRichardson, T. & Fairbairn, M. 2014, MNRAS, 441, 1584  \nRidolfi, A., Gautam, T., Freire, P. C. C., et al. 2021, MNRAS, 504, 1407  \nShklovskii, I. 1970, Soviet Astronomy, Vol. 13, p. 562, 13, 562  \nSkilling, 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, & U. V. Toussaint, 395-405  \nSkilling, J. 2006, Bayesian Analysis, 1, 833  \n- Smith, P. J., Hénault-Brunet, V., Dickson, N., Gieles, M., & Baumgardt, H. 2024, arXiv e-prints, arXiv:2407.06274\n- Sollima, A., Bellazzini, M., Smart, R. L., et al. 2009, MNRAS, 396, 2183\n- Speagle, J. S. 2020, MNRAS, 493, 3132\n- Strigari, L. E., Bullock, J. S., & Kaplinghat, M. 2007, ApJ, 657, L1\n- Trager, S. C., King, I. R., & Djorgovski, S. 1995, AJ, 109, 218\n- van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., & de Zeeuw, P. T. 2006, A&A, 445, 513\n- van der Marel, R. P. 1994, MNRAS, 270, 271  \nVan der Marel, R. P. & Anderson, J. 2010, The Astrophysical Journal, 710, 1063  \n- Vasiliev, E. & Baumgardt, H. 2021, MNRAS, 505, 5978  \nVerbunt, F. & Freire, P. C. C. 2014, A&A, 561, A11  \nVirtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261  \n- Vleeschower, L., Corongiu, A., Stappers, B. W., et al. 2024, MNRAS[ arXiv:2403.12137 ]  \nWang, L.-C. & Xie, Y. 2021, Research in Astronomy and Astrophysics, 21, 270  \n- Watkins, L. L., van de Ven, G., den Brok, M., & van den Bosch, R. C. E. 2013, MNRAS, 436, 2598\n- Watkins, L. L., van der Marel, R. P., Bellini, A., & Anderson, J. 2015a, ApJ, 803, 29\n- Watkins, L. L., van der Marel, R. P., Bellini, A., & Anderson, J. 2015b, ApJ, 812, 149\n- Wilkinson, M. I., Kleyna, J., Evans, N. W., & Gilmore, G. 2002, Mon. Not. Roy. Astron. Soc., 330, 778\n- Wirth, H., Bekki, K., & Hayashi, K. 2020, Monthly Notices of the Royal Astronomical Society: Letters, 496, L70  \nXie, Y. & Wang, L.-C. 2020, Research in Astronomy and Astrophysics, 20, 191 Xie, Y., Yin, D., Wang, L., et al. 2024, MNRAS, 527, 7743  \nYin, D., Zhang, L.-y., Qian, L., et al. 2024, arXiv e-prints, arXiv:2405.18228  \nZentner, A., Dandavate, S., Slone, O., & Lisanti, M. 2022, JCAP, 07, 031  \nZhang, L., Freire, P. C. C., Ridolfi, A., et al. 2023, ApJS, 269, 56  \nZhao, H. 1996, MNRAS, 278, 488  \nZhao, J. & Heinke, C. O. 2022, MNRAS, 511, 5964  \nZhao, J. & Heinke, C. O. 2023, MNRAS, 526, 2736  \nZocchi, A., Gieles, M., & Hénault-Brunet, V. 2017, MNRAS, 468, 4429  \nZocchi, A., Gieles, M., & 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 software packages: dynesty (Speagle 2020; Koposov et al. 2022), GravSphere (Read & Steger 2017; Read et al. 2018; Genina et al. 2020; Collins et al. 2021), corner.py (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': '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  \n<!-- image -->  \n∞  \nFig. 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 D: Use of different dark components': '<!-- image -->  \n/circledot  \n/circledot  \nFig. 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.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.  \n<!-- image -->  \nFig. 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.  \n<!-- image -->  \n† 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.  \ncomponent 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).  \nFig. D.2. Posterior distributions for the mass model parameters using a Gaussian profile (Eq. (D.1)) for the extended dark component.  \n<!-- image -->  \n<!-- image -->  \nFig. C.2. Combined velocity dispersion profiles for the three components using the Anderson & Van der Marel (2010) center (And, left) Noyola et al. (2010) (Noy10, middle) and Noyola et al. (2008) (Noy08, right) kinematic centers. The maximum-posterior velocity dispersion profile (red) from our main analysis using the And center is shown for reference in all plots, yielding fits of comparable quality for all the cases considered. Due to the close-to-isotropic behavior of our inferred distribution at the 5.2 kpc distance employed, the three components of our maximum posterior fits are almost identical, with the large majority of the data points showing close overlap with each other.  \n<!-- image -->  \nFig. D.1. Posterior distributions for the mass model parameters using the generalized mass profile in Eq. (17) for the extended dark component.  \n<!-- image -->  \nFor the second case we adopt the density of a spherical Gaussian profile as in Evans et al. (2022), which can be expressed as:  \nρ cen , g( r ) = M cen (2 π r 2 g ) 3 / 2 exp -r 2 2 r 2 g ! . (D.1)  \nIn 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 kinematically relevant IMBH component, capped at a few thousand solar masses.  \n<!-- image -->  \n/circledot'}