Timing of millisecond pulsars in NGC 6752 - III. On the presence of non-luminous matter in the cluster’s core.

Millisecond pulsars are subject to accelerations in globular clusters that manifest themselves in both the ﬁrst and second spin period time derivatives, and can be used to explore the mass distribution of the potentials they inhabit. Here we report on over 20 years of pulsar timing observations of ﬁve millisecond radio pulsars in the core of the core-collapse globular cluster NGC 6752 with the Parkes (Murriyang) and MeerKAT radio telescopes, that have allowed us to measure the proper motions, positions and ﬁrst and second time derivatives of the pulsars. The pulsar timing parameters indicate that all the pulsars in the core experience accelerations and jerks that can be explained only if an amount of non-luminous mass of at least 2 . 56 × 10 3 M ⊙ is present in the core of NGC 6752. On the other hand, our studies highly disfavour the presence of an intermediate mass black hole at the center of the cluster, with a mass equal to or greater than ∼ 3000 M ⊙ .


Introduction
Globular clusters (GCs) are among the oldest and densest structures in a galaxy and often contain several ×10 5 stars in spherical structures with radii of just a few tens of parsecs.In recent decades, more and more powerful optical telescopes have been investigating these systems, providing new insights into their mass distribution and internal structure (see, e.g., Baumgardt & Hilker 2018;Leitinger et al. 2023, and references therein).Today accurate censuses of their stellar populations can be obtained via optical data, and their structure, dynamical status, and evolution can be explored with N-body simulations.While their macroscopic properties are increasingly well understood, their dense cores remain difficult to probe.Even high spatial resolution instruments like GAIA (Gaia Collaboration et al. 2016) struggle to resolve stars in the central regions due to the high number densities and distances.Another reason is the presence in the central regions of high-mass nonluminous stellar remnants like neutron stars (NS) and possibly stellar mass black holes (Bahcall & Wolf 1976, Bahcall & Wolf 1977), and also massive underluminous white dwarfs (e.g., Vitral et al. 2023) which have migrated to the clusters' cores, thus leaving the lightest objects at larger distances from the center due to mass segregation.This scenario finds confirmation via the observed spatial distribution of millisecond pulsars (MSP) in GCs, which are on average among the heaviest objects in a cluster and usually are located near the cores.The pulsars located outside the core have been postulated to be the result of multibody ejection mechanisms, some involving the presence of an intermediate-mass black hole (IMBH) of several tens of solar masses (e.g., Colpi et al. 2002), located near the cluster's center of gravity.In recent works, instead, it has been argued that three-body interactions involving stellar mass objects are a viable mechanism for expelling binary MSPs from the cluster core (see, e.g., Leigh et al. 2024).
The MSPs in the GC NGC 6752 show characteristics that are typical of the populations of core-collapsed GCs: a large predominance of isolated pulsars (for other examples, see the pulsar population of M15, Anderson 1993, NGC 6624, Abbate et al. 2022and NGC 6517, Pan et al. 2021) and a much larger fraction of pulsars at large distances from the core.As discussed in detail by Verbunt & Freire (2014), these facts are likely causally related: the large number of stellar interactions per binary in the core region-to which NSs migrate via dynamical friction-implies that many MSP binaries are disrupted in three-body encounters, with many being ejected from the cluster, and with several others barely staying bound to it.The ejections of NSs are more likely if they are involved in interactions with more massive stellar remnants.
However, even in comparison with the pulsar populations of other core-collapsed GCs, the population of NGC 6752 is extreme in these respects.Its peculiarities are hardly explainable without invoking the presence of a nonnegligible amount of under-or nonluminous (hereafter, we name all of that as nonluminous) mass in the core, under the form of massive remnants.PSR J1910-5959A (hereafter PSRA, D' Amico et al. 2001), the first pulsar discovered in this cluster and its only binary pulsar known to date, has the largest offset from the center of its host GC in core radius units, namely ∼58 r c (this fact has cast doubts on the association, which has, however, been demonstrated recently by Corongiu et al. 2023) and angular units (~¢ 6.3); while PSR J1910−5959C (hereafter PSRC, D 'Amico et al. 2002, hereafter Paper I) has the third angular offset in core radius units, namely ∼23 r c , and fourth in angular units (∼2 69 ).For these two pulsars, recent episodes of ejection from the core are a very likely possibility, since dynamical friction should have taken them back to the core in less than ∼1 Gyr (Colpi et al. 2002).The other known pulsars, namely PSR J1910-5959B, PSR J1910-5959D, and PSR J1910-5959E (hereafter PSRB, PSRD, and PSRE respectively, Paper I), and the recently discovered PSR J1910−5959F (hereafter PSRF, Ridolfi et al. 2021) are located in the core of NGC 6752.We will refer to these four pulsars as the core pulsars.
In Paper I it was shown that the measured spin period derivatives of PSRB, PSRD, and PSRE are dominated by the acceleration imparted to them by the cluster potential well, and at least 1.3 × 10 4 M e should be present in the form of nonluminous objects in the core.The measured spin period derivative of PSRF (Ridolfi et al. 2021) is of the same order of magnitude as the other three core pulsars, and also its value is likely due to the gravitational effect of nonluminous matter in the core of NGC 6752.
In this work, we report updated timing solutions for PSRB, PSRC, PSRD, PSRE, and PSRF10 based on observations taken over ∼20 yr with the Parkes Murriyang and MeerKAT radio telescopes, and investigate the total mass and distribution of the nonluminous matter in the cluster core.This paper is organized as follows: In Section 2 we describe the observations and the data analysis, and in Section 3 we comment on the proper motions of the core pulsars.In Section 4 we place constraints on the nonluminous matter in the cluster core, and in Section 5 we discuss our results in terms of the central mass-to-light ratio.
In Section 6 we comment on the main sources of uncertainty in our analysis, and we review our numerical results.Finally, in Section 7 we summarize our findings.

Parkes Observations
The GC NGC 6752 was observed with the Parkes Murriyang 64 m radio telescope, located in Australia, from 1999 September to 2016 March.Observations were carried out in the L-band with the Multibeam and the H−OH receivers, depending on availability, with a bandwidth of 256 MHz and using two orthogonal polarizations that were detected and summed before digitization.The Stokes I signal was processed with the 1 bit Analog Filterbank (AFB; see, e.g., D 'Amico et al. 2001) digital signal processor (hereafter "backend") up to its decommissioning in 2012 September.The AFB acquired the detected signal at the central frequency of 1390 MHz, split the observed bandwidth into 512 frequency channels, and recorded the data with a sampling time of 80 μs (125 μs in some earlier observations).In all observations taken after 2012 September, and in some before, the data were recorded in search mode with the Pulsar Digital Filterbank 3 (PDFB3) and 4 (PDFB4; see, e.g., Manchester et al. 2013 and references therein).These backends acquired the detected signal at the central frequency of 1369 MHz, split the observed bandwidth in 512 (DFB3) or 1024 (DFB4) frequency channels, and recorded the data with a sampling time of 80 μs (125 μs in a few observations).
We reduced and analyzed the raw data with the same method used and explained in detail in Paper I and Corongiu et al. (2006, hereafter Paper II).Briefly, for each pulsar, we folded the data at the approximate topocentric period with the software package dspsr11 (van Straten & Bailes 2011), as predicted by the ephemeris published in Paper II.We folded the data with a typical subintegration length of 1 minute and kept the same number of frequency channels as the raw data to aid in accurate dedispersion.The resulting reduced data, whose files are usually described as "timing archives," have the form of a matrix whose elements are the binned profile amplitudes observed in each frequency channel, one for each subintegration time, and some metadata (header) information.

MeerKAT Observations
The MeerKAT radio telescope (Jonas & MeerKAT Team 2016;Camilo et al. 2018), located in South Africa, is an array composed of 64 offset-Gregorian dishes, each equivalent to a circular aperture of approximately 13.5 m across.We used MeerKAT and its UHF-band (544-1088 MHz) and L-band (856-1712 MHz) receivers to observe NGC 6752 on several occasions between 2019 July and 2022 January, as part of the MeerTime12 (Bailes et al. 2020) and the TRansients And PUlsars with MeerKAT (TRAPUM13 ; Stappers & Kramer 2016) Large Survey Projects.Each of these projects has its own backend, tailored to its main scientific goals.MeerTime, which is mostly devoted to precision pulsar timing and polarimetry, makes use of the Pulsar Timing User Supplied Equipment (PTUSE) backend.This is capable of acquiring up to four Nyquist-sampled tied-array beams with up to 4096 coherently dedispersed frequency channels with four Stokes parameters.The data were recorded as 8 bit search mode PSRFITS files with a time resolution of 9.57 μs at L-Band and 7.53 μs at UHF.Given that only pulsars A and C are located far from the core of the cluster, while all the others are less than 0 2 away from the nominal center of the cluster, three PTUSE beams (centered on pulsar A, on pulsar C, and on the nominal center of NGC 6752, respectively) were sufficient to record the signals from all the known pulsars.Also, being the DMs of the pulsars all in the range of 33.20-33.70pc cm −3 , a single search mode file coherently dedispersed at a DM = 33.3pc cm −3 was sufficient to allow the redetection of all the pulsars without any loss of sensitivity.TRAPUM, on the other hand, is focused on the search for pulsars and radio transients, and uses two computing clusters.First, the Filterbanking BeamFormer User Supplied Equipment (FBFUSE) cluster applies the beamforming technique to combine the raw signals from all the antennas and synthesize hundreds (up to 288 for our NGC 6752 observations) tied-array beams on the sky (Barr 2018;Chen et al. 2021).These are then recorded as search mode "filterbank" files by the Accelerated Pulsar Search User Supplied Equipment (APSUSE) cluster: the observing band is recorded with a typical time resolution of ∼60-80 μs and is split into 4096 frequency channels.The fine channelization is particularly important in this case to remove the effects of interstellar dispersion, as APSUSE does not benefit from coherent dedispersion, due to computational constraints.The majority of the MeerKAT observations of NGC 6752 were carried out using PTUSE and APSUSE simultaneously, so as to take advantage of their complementary characteristics.The exact used setups depended on the main scientific purpose of each observing session, as well as on the availability of new observing modes (such as the possibility of recording more than a single beam with PTUSE, or to record a different number of channels between PTUSE and APSUSE during a simultaneous MeerTime+TRAPUM session) which were gradually implemented over the course of these projects.We excluded the APSUSE data for PSRB, because the time stamp was not correctly recorded in the corresponding data files.We report all the MeerKAT observations of NGC 6752 used in this paper, along with their configurations, in Table 1.

ToA Extraction and Timing
We extracted the pulse times of arrival (ToAs) with the routines of the software suite psrchive (Hotan et al. 2004), by coherently adding in phase the profiles in each archive with respect to subintegrations and channels, and convolving them with a high signal-to-noise ratio (S/N) template obtained by summing in phase the observed profiles with the best S/N.For each pulsar, we produced a specific template for each telescope-backend-observing mode combination.The pulsars were often subject to strong scintillation, due to the ionized interstellar medium, which produces random variations of the pulse S/N with respect to both time and frequency, correlated on some timescale.For this cluster the timescale is such that amplitude variations occur not only from one observation to another, but also within single observations, and they are so dramatic that pulses can be very bright at some moments and/ or in some frequency subbands, and completely undetectable in others.For these reasons, we could not determine a single optimal integration length and subband frequency width along which the profiles in the archive should always be summed; instead, we visually inspected every single archive to determine the best time interval(s) and frequency band(s) where pulses were evaluated as detections.Whenever possible, we extracted more than one ToA from a single observation, either by decimating in time, frequency, or both.
We determined the instrumental time jump between the different backends used for Parkes observations, and between the Parkes and MeerKAT data sets as a whole.We determined these jumps for each pulsar separately, since they are not due to technical reasons only, but mainly because templates differ, from one pulsar to another, in shape, displacement, and time extent.We used the ToAs obtained from the Parkes AFB backend as reference (i.e., we determined the time jumps with respect to), since for all pulsars these sets span the largest time interval and contain the largest number of ToAs.At first we considered the ToAs for each telescope separately, and fitted them against the best timing model for each set.We took the resulting χ 2 and multiplied the ToA uncertainties by its square root.Once this operation is performed, a new fit on the same ToAs against the same timing model should return a reduced χ 2 = 1.We then combined the two telescopeʼs ToAs and fitted them for a single time jump alone, using Parkes ToAs as the reference.Finally, we proceeded with the fit against the pulsar timing model, still allowing this time jump to vary for correctly estimating the uncertainties in the timing modelʼs parameters.
We fitted the obtained ToAs against the rotational and astrometric parameters with the pulsar timing software TEMPO2 (Hobbs et al. 2006), using the DE430 solar system ephemeris from the NASA Jet Propulsion Laboratory14 for barycentering.Table 2 reports our timing results.The values for each parameter in the timing model are reported with the nominal TEMPO2 1σ uncertainty on the last significant digit (in parentheses).The main improvement in the timing solutions with respect to Paper II is the determination of the second derivative of the spin period for all pulsars, the proper motions for PSRB, D, and E, and the refinement of the latter for PSRF.Parameters whose measurement has already been reported in previous works are now determined with a similar or higher precision.Figure 1 displays the corresponding timing residuals.The amplitude of the plotted error bars is equal to the measured one for each ToA, multiplied by the amount necessary to obtain χ 2 = 1 in the preliminary fit of each telescope data set separately (see above).No clear trend is evident in the   residuals of any pulsar, thus indicating that the adopted timing models adequately describe the available data sets.

Proper Motions of the Core Pulsars and the Central Velocity Dispersion
The pulsars in the core of NGC 6752 have excellent proper motion determinations, namely with a significance of 11σ, 25σ, 31σ, and 11σ for PSRB, PSRD, PSRE, and PSRF in alphabetical order, respectively.The average motion of the pulsars results in m d á ñ = - a cos 3.15 0.51 mas yr −1 (quoted uncertainties are everywhere at the 1σ level unless explicitly indicated) and 〈μ δ 〉 = −4.13± 0.31 mas yr −1 , in excellent agreement with the cluster motion m d = - a cos 3.162 0.023 mas yr −1 and μ δ = −4.028± 0.022 mas yr −1 , obtained from GAIA optical observations (all clusterʼs parameters picked up from literature and used in this work are presented in Table 3 with all necessary references, unless explicitly indicated).At the cluster distance of 4.125 ± 0.040 kpc the relative projected transverse velocity of each pulsar with respect to the cluster center, for which the GAIA proper motion is assumed, are 15.5 ± 4.7, 7.7 ± 4.1, 4.8 ± 3.6, and 15.4 ± 11.1 km s −1 respectively, thus implying that (barring a large radial velocity component) all pulsars are compatible with being bound to the cluster core, since the escape velocity from the core is V esc = 34.5 km s −1 .
We estimated the 1D velocity dispersion of the core pulsars, by calculating the standard deviation of the proper motion components in R.A. and decl.separately, then averaging the values in the two directions.We obtain σ V,PSRs = 8.0 ± 1.9 km s −1 .This value is compatible with the value σ 0 = 8.5 ± 0.2 km s −1 , resulting from the measured proper motion dispersion in the core of the cluster, σ μ = 0.436 ± 0.009 mas yr −1 , using optical observations.Two pulsars in the core of the cluster, PSRB and PSRE, have a negative measured time derivative of their spin period  P ; meas that is the signature of their accelerated motion in the gravitational potential well of the cluster (Paper I), since  P is always intrinsically positive because of pulsar spindown.In the case of PSRD (Paper I) and PSRF (Ridolfi et al. 2021, this work)  P meas is positive, with an absolute value at least as large as for PSRB.An insight into the Galactic population of MSPs15 (here defined as pulsars with a spin period not larger than 10 ms) shows that it is highly unlikely that  P meas for PSRD and PSRF is dominated by their intrinsic spindown  P intr .In fact, the Galactic population of MSPs consists of 228 objects for which it has been measured  > P 0. Among them, 223 MSPs have  < - P 10 19 s s −1 , with an average spindown rate  á ñ = ´-P 1.40 10 20 s s −1 , and a standard deviation ´-1.22 10 P 20 s s −1 .For the four pulsars in the core,  P meas | |is larger than  á ñ P by factors that range from 31 (in the case of PSRE) to 69 (in the case of PSRD).The maximum  P in this sample is 7.74 × 10 −20 s s −1 (PSR J0218+4232, Desvignes et al. 2016), a value that is an order of magnitude smaller than  P meas of PSRF, namely the smallest value among the core pulsars for which  > P 0 meas .For four of the remaining five Galactic MSPs,  P is comprised between 1.05 × 10 −19 s s −1 (PSR B1937+21, Reardon et al. 2021) and 1.63 × 10 −19 s s −1 (PSR J1850+0242, Scholz et al. 2015).But even in the case of PSR J1850+0242, the  P value is still a factor of ∼4.5 smaller than  P meas for PSRF.Moreover, if one includes these four MSPs in the aforementioned calculation of the average spindown, one obtains  á ñ = ´-P 1.60 10 20 s s −1 and  s =

P
´-2.00 10 20 s s −1 .It is immediately seen that the  P of these four MSPs is more than 4σ away from the average.The fifth one is PSR J1402+13 (P = 5.9 ms,  = ´-P 4.8 10 17 s s −1 , Abdollahi et al. 2022), about which no discussion seems to be present in the literature.Because these five high  P objects represent the 2.2% only of the sample of Galactic MSPs, and their  P is significantly higher than the average MSPs spindown rate, we evaluated these five objects as outliers in the  P distribution of the Galactic MSPs, and assumed that the distribution of the intrinsic spindown of the galactic MSPs is well represented by the 223 objects for which  < < - P 0.0 10 19 s s −1 .It could be objected that an MSP in a GC might have a  P intr that is much higher than the average of the Galactic MSPs, since the clusterʼs crowded environment might have induced a different formation path for the hosted MSPs, which results in a considerably large  P intr .Nevertheless, NGC 6752 also hosts PSRA, for which Corongiu et al. (2023) obtained  = ´-P 5.02 10 intr 21 s s −1 , and PSRC, whose  P meas and projected angular separation from the clusterʼs center of gravity imply   ´-P 1.71 10 intr 20 s s −1 .Since for both PSRA and PSRC  P intr is fully consistent with the average of the Galactic field MSPs, we can safely assume that this is also true for the pulsars in the core of NGC 6752.
The effects of the cluster environment on the core pulsars are also evident in the second time derivative of their spin periods.
In fact, the upper limits, both observational and theoretical, on the intrinsic P ̈for MSPs are so low, that it can safely be Central density (M e pc −3 ) Central proper motion dispersion (mas yr −1 ) assumed that the intrinsic P ̈for the core pulsars are completely negligible in comparison to the measured ones, whose origin can thus be attributed to the clusterʼs dynamics.
In Paper I it was investigated whether the acceleration imparted to the core pulsars (known at that time) could be entirely due to the observed luminous mass, or if a further amount of matter was necessary to explain their  P meas .The authors of Paper I used the argument of the mass-to-light ratio    in the central regions (see Equation (1) in Paper I), and found that a minimum    of an order of 10, namely >    9 for PSRB and E, and 13 for PSRD, is required to impart to these pulsars the accelerations, whose component along the line of sight produces the measured  P. Hereafter in this paper, in order to improve the reading of the text, the term acceleration will indicate the component along the line of sight of an acceleration, unless explicitly specified otherwise.These    values were much larger than the one available at that time and obtained with optical observations, =    1.1 (Pryor & Meylan 1993).The authors of Paper I calculated the amount of mass that should be present in the core of the cluster in the form of nonluminous objects, and found a lower limit of 1.3 × 10 4 M e .

Our Approach
Given the high precision of our determination of both the first and second derivatives of the spin period for the core pulsars, we revisited the estimate of the nonluminous mass in the core of NGC 6752 with a different approach.Phinney (1993) detailed how  P meas is related to the intrinsic spindown  P intr and all accelerations acting on a pulsar in a GC with the following relation 16 : where c is the speed of light, a l,MW is the acceleration imparted by the Milky Way,17 a l,SHK is the apparent acceleration due to the pulsar proper motion, also known as the Shklovskii (1970) effect, and a l,GC is the acceleration imparted by the GC.When the second derivative of the spin period is also measured, a further equation is available, that links the measured P ̈to the time derivative of the imparted acceleration, often referred to as the jerk.In fact, taking the time derivative of Equation (1) and considering relevant terms only, one obtains:

̈̈( )
where we have introduced the effect of the nearest neighbor star,  a l,NN , which cannot be in general neglected in crowded environments like GCs (see, e.g., Abbate et al. 2019), while we dropped the time derivatives of the Milky Way acceleration and that of the Shklowskii effect.These two terms can be neglected because, across the epoch range spanned by the observations, the 3D motion of the cluster in the Milky Way does not produce significant changes on its position, hence on the acceleration imparted by the Galaxy.The pulsar proper motion changes even less, hence the apparent acceleration due to its transverse motion remains substantially constant.
The acceleration and jerk imparted by the cluster can thus be determined by solving Equation (1) with respect to a l,GC and Equation (2) with respect to  a l,GC : and the link with the cluster structural parameters is obtained by replacing a l,GC in Equation (3) and  a l,GC in Equation (4) with the expression that results by assuming a given mass distribution for the cluster.In the simplest case of a stationary and spherically symmetric distribution for the mass of the cluster, a l,GC takes the well-known simple form: where M(r) is the cluster mass enclosed in a sphere of radius r, centered at the center of gravity of the cluster, and G is the Gravitational constant (6.67 × 10 −11 m 3 kg −1 s −2 ).We express a generic position, with respect to the cluster center, in terms of two coordinates perpendicular to the line of sight, a a d º - 2 obviously holds.In these definitions, α PSR and δ PSR are the celestial coordinates of a given pulsar, while α GC and δ GC are the celestial coordinates of the cluster gravity center, and r l is defined so that its value is positive for positions farther than the cluster center with respect to the observer, and negative otherwise.At this stage, we considered r α and r δ of each pulsar as known exactly, despite the uncertainties on the assumed coordinates of the center of gravity of the cluster (0 5 in both coordinates, Ferraro et al. 2003).In Section 6 we will discuss the impact of these uncertainties on our results.
The analytic expression for  a l,GC is obtained by taking the time derivative of Equation (5): where v is the pulsar velocity with respect to the cluster center, whose components are GC and μ δ,GC are the proper motion components of the cluster, and v l ≡ dr l /dt.Given the above definition of r l , v l is positive when the pulsar is moving away from the observer, and negative otherwise.We did not consider the quantities v α and v δ to be exactly known, since their uncertainty is comparable to their amplitude and to the central velocity dispersion.
The problem is so far underdetermined: for a given pulsar, only two equations are available, but the unknown terms are its depth inside the cluster r l , its 3D velocity v, and all the parameters that enter in the analytic expression of a given mass distribution for the cluster.We indicate the set of all of these parameters with the formal vector s ≡ (s 1 , s 2 ,K,s N ).Moreover, the contribution to  a l,GC due to the nearest neighbor calls into play some further parameters that describe the cluster structure and dynamic properties at the pulsar position (see below).

The Bayesian Analysis
The availability of more than one pulsar does not allow us to fully constrain the problem, since each object has its own position and velocity, and also, the clusterʼs local properties differ from a given position to another.For this reason we adopted a Bayesian approach.All terms in the right-hand side of Equation (3) but the intrinsic spindown (see below) are known.We can then assume that they obey a Gaussian probability distribution with a mean value equal to their value, either measured or derived, and a standard deviation equal to their 1σ uncertainty.This in turn implies that also a l,GC (hereafter a l for notation simplicity) follows a Gaussian probability distribution: whose mean value a l,m is given by Equation (3), and whose standard deviation s a l,m is given by the usual rules for the propagation of the uncertainties.Once a l is substituted with the explicit expression given by Equation (5), where in turn M(r) is substituted with the analytical expression that corresponds to a chosen model for the mass distribution, Equation (7) returns the probability distribution ), for the structural parameters in s and the still unknown pulsar depth in the cluster r l , as it results by considering a single pulsar alone, identified by the formal index X (the subscript a l in the symbol s P r , a l ,X l ( ) indicates that this probability distribution comes from the analysis of the acceleration, in order to distinguish it from the distribution that results from the analysis of the jerk).
We calculated the Milky Way acceleration by applying Equation (16) in Lazaridis et al. (2009), but using for the vertical component of the Galactic acceleration F z the analytic formula provided by Li & Widrow (2021) (Equation ( 14 )to the uncertainty on a l,m that results from the standard propagation of the uncertainties on the three known terms.
We applied the method illustrated above to Equation (4), for deriving the probability distribution obeyed by  a l,GC (hereafter  a l for notation simplicity), which in turn translates into the probability distribution  P a ,X l for the structural parameters in s and the unknown quantities for a generic pulsar X.In this case it must be taken into account that the nearest neighbor contribution on the measured jerk cannot be calculated but only statistically treated and that, most important, it does not follow a Gaussian distribution, but a Lorentzian one (see, e.g., Abbate et al. 2019 and references therein): The quantity  a 0 is dubbed as characteristic jerk, and it is related to the clusterʼs structure (Prager et al. 2017) by the following relation: where ξ ; 3.04 is a numerical constant, while σ v , n, and á ñ m are the clusterʼs 1D velocity dispersion, star number density and mean stellar mass at the pulsar position.The product á ñ n m can be set equal to the mass density at the pulsar position, while σ v can be only constrained to be smaller than the central 1D velocity dispersion.Once a model is assumed for the clusterʼs mass distribution, the probability distribution for  a NN depends on the pulsar position and the local velocity dispersion σ v .Therefore,  P a ,X l depends on the parameters in s, the pulsar depth in the cluster, the pulsar 3D velocity in the cluster frame, and the clusterʼs 1D velocity dispersion at the pulsar position.

Because of the peculiar distribution for
) must be calculated from more general principles of statistics.Let α and β be two generic quantities, which follow probability distributions P α (α) and P β (β) respectively, and γ = α − β their difference.The probability distribution for γ, P(γ), is given by the following integral 18 : In our case a = cP P ̈, β is indeed  a NN , and γ is the value for  a l as predicted by Equation (6) once a given mass distribution is assumed.Because of the measurement of the second-order time derivative of the spin period, P ̈follows a Gaussian distribution with mean value P meas ̈, and standard deviation s P meas ̈.As a direct consequence, once P is kept fixed, also cP P ̈follows a Gaussian distribution with mean value cP P meas ̈and standard deviation s c P P meas ̈.We completely neglected the contribution due to the intrinsic evolution of the spin period, namely P intr ̈since, as we already commented on, it is negligible in comparison with P meas ̈.The joint probability of simultaneously having an acceleration a l and a jerk  a l acting on a given pulsar X is the product of the two above-mentioned distributions.Its marginalization with respect to r l , v and σ v returns the probability distribution P X (s) for the parameters that describe the assumed mass model: 18 One can easily demonstrate that, if α and β follow Gaussian distributions centered at α 0 and β 0 with standard deviations σ α and σ β respectively, Equation (10) results in another Gaussian function centered on α 0 − β 0 , and with standard deviation s s + a b 2 2 , i.e., the usual result of the uncertainty propagation for the difference of two Gaussian distributed quantities.
where the P ˜functions are the prior probability distributions for their arguments.We assumed that P r l ˜( ) is flat and nonzero only for those values for r l that place a given pulsar within one core radius from the cluster center, i.e.,   r r r 0 , and whose sign is opposite to  P meas , namely r l > 0 if  < P 0 meas and vice versa.In fact, both a l,MW and a l,SKH are positive, and their sum is so smaller than  c P P meas | | for all core pulsars, that Equation (1) is satisfied only if a l,GC has the same sign of  P meas , and the explicit expression for a l,GC , as given by Equation (5), states that the sign of r l must be opposite to a l,GC , hence to  P meas .We assumed v P ˜( ) to be the product of the priors a P v ˜( ), d P v ˜( ), and P v l ˜( ) for the three components of the velocity.Thanks to our measurement of the proper motion for all core pulsars, we assumed that a P v ˜( ) and d P v ˜( ) are Gaussian distributions, whose mean and standard deviation can be obtained from their definition and the usual uncertainty propagation rules.Because no information is available about the radial velocity of the core pulsars, we assumed P v l ˜( ) to be flat and nonzero in the range −50 km s −1 v l 50 km s −1 , where the limits are a conservative rounding of the escape velocity from the core, and we allowed both positive and negative values since it is not a priori known whether a pulsar is moving away or toward the observer in the cluster frame.Finally, we assumed s P v ˜( ) to be flat in the range 0 σ v σ 0 , since it cannot be larger than the central velocity dispersion.Its upper limit will be discussed later, since it is model dependent.
Because the same mass distribution must be responsible for the accelerations and the jerks observed in all the core pulsars, the final probability distribution P(s) is given by the product of distributions P X (s), later multiplied by the prior distribution s P ˜( ):

The Mass Models and the Distribution of the Luminous Mass
The models that we considered for the mass distribution in the core of the cluster are the sum of two components.The first one is due to the luminous mass, whose spatial distribution is described by the model that best fits to the data taken with optical observations.We assumed it was relatively well known and thus we kept it fixed.The second component is due to the nonluminous mass, and we investigated two scenarios for its distribution.In the first one we simply assumed that all the nonluminous mass is generically contained in a sphere, whose radius is smaller than the true distance of the closest pulsar to the cluster center (see Section 4.4).This scenario allowed us to obtain a first estimate of the amount of nonluminous matter in the core of the cluster.In the second scenario, instead, we assumed that an IMBH is located at the center of gravity of the cluster.This scenario differs from the first one (see Section 4.5), since the presence of an IMBH has a substantial impact on the mass distribution in the core, and it allowed us to test whether the resulting picture represents a physically meaningful description of the mass distribution in the core of the cluster.
The distribution of the luminous mass obviously plays a key role in our work.In fact, if one adopts a model that, e.g., underestimates its content, one obtains an overestimation of the nonluminous mass, and vice versa.For this reason, we paid great attention to the previous works, available in the literature, on the mass distribution of NGC 6752, and we identified a model that allows one to obtain a consistent picture of the optically observed structure of this cluster.
The most recent study of the mass surface and radial profile of NGC 6752 was reported in Baumgardt & Hilker (2018).The authors performed N-body simulations of the formation and evolution of the Galactic GCs.The results of their simulations for NGC 6752 (see their Figure E13) are consistent with the data at angular radii θ ⊥  10″, but they clearly underestimate the observed surface density profile at smaller radii.This evidence can be considered consistent with the results obtained by Ferraro et al. (2003), who studied the surface star density profile of NGC 6752 using different models.In particular, they obtained a very good description of their data at all angular radii by using a combination of two King (1966) density profiles 19 (see text and Figure 5 in Ferraro et al. 2003), for separately describing the cluster star density in its inner and outer regions, in projection on the sky.We recall that a King profile describes the structure of a GC in terms of its mass density, whose analytic expression is given by the following equation: where ρ 0 is the central mass density and r c is the core radius.It is interesting to note that the transition between the two profiles occurs at θ ⊥ ∼ 10″, i.e., at about the same angular radius at which the simulations by Baumgardt & Hilker (2018) begin to agree with the data.Since we have assumed that all core pulsars are indeed located in the core of the cluster, i.e., their true 3D distance from the center of gravity of the cluster is smaller than the core radius, which in turn is smaller than the angular radius at which the outer King profile begins to describe the data, we can safely assume that the luminous mass distribution that contributes to the acceleration imparted on the pulsars is well described by the inner King profile only.We tested whether the the double King profile by Ferraro et al. (2003) adequately describes the cluster profile at all radii, by calculating the cluster total mass that results from all the above-mentioned assumptions, and the already known clusterʼs parameters (see Table 3).We set the value for the core radius of the outer King profile to r c = 28″, according to the results by Ferraro et al. (2003), and we normalized it by imposing that at r = 10″ it must return the same density given by the inner one.Finally, we integrated this two-component density profile up to the radius r NGC 6752 = 12.26 pc, which corresponds to 4.27 projected half-light radii, thus following the prescription given by Leitinger et al. (2023; see their Section 4.1 for a discussion about this choice).We obtained a total mass of the cluster M TOT = 2.865 × 10 5 M e , a value just a few percent higher than the value (2.76 ± 0.04) × 10 5 M e (see Table 3), thus ensuring that the double King profile adequately describes the optically observed mass profile of the cluster at all radii.

A First Estimate of the Nonluminous Mass in the Core
Our first estimate of the nonluminous mass M NL in the cluster core is based on assuming that it is entirely contained in a sphere whose radius is smaller than the distance of the closest pulsar to the cluster center.We simplified the problem by assuming that the nonluminous mass does not perturb the distribution of the luminous matter.With these assumptions, this distribution for the nonluminous mass can be mathematically treated as a point-mass M NL placed at the center of gravity of the cluster.For this reason, this model will later be referred to as the point-mass model.The formal vector s is thus the 1D vector s ≡ (M NL ).We assumed the prior distribution P M NL ˜( ) to be flat and nonzero in the range 0 M NL 10 4 M e .In this model we set to 10 km s −1 the upper limit for the local velocity dispersion that enters in the characteristic jerk expression.This value is a conservative upper limit for the measured central velocity dispersion, namely σ 0 = 8.5 ± 0.2 km s −1 .We obtained at the 2σ and 3σ level, respectively).This means that about 3 thousand solar masses of nonluminous mass, with a 3σ lower limit M NL 2.68 × 10 3 M e , are required to explain the measured accelerations and jerks acting on the core pulsars.

Is There an IMBH in the Core of NGC 6752?
We explored a scenario where an IMBH resides at the center of gravity of NGC 6752.The presence of an IMBH modifies the density profile of a GC, with respect to the predictions of a pure King profile, up to a radius r i dubbed as the IMBH influence radius (Baumgardt et al. 2004a and references therein).Such a density distribution is called cusp, since for r r i it predicts higher values for the density than the ones predicted at the same radii by the King model that describes the density profile at r > r i .In this scenario the nonluminous mass is not only due to the IMBH alone, but also to the extra mass M CUSP − M King (r i ), namely the difference between the amounts of the cluster mass contained up to r = r i , as predicted by the cusp density law (M CUSP ), and the King model that describes the luminous mass at radii r r i (M King (r i )).Prager et al. (2017) explicitly obtained the analytic expression for the density profile of the cusp ρ C (r), as it results from the works by Baumgardt et al. (2004aBaumgardt et al. ( , 2004b)): The influence radius r i is related to the IMBH mass M BH and the 1D velocity dispersion of the stars in the cluster core σ 0 by the formula (Baumgardt et al. 2004a, Equation (1)): Outside r i the mass density profile follows the standard King (1966) law (see Equation ( 13)).In the simulations by Baumgardt et al. (2004b) the radial density profile is a continuous function at r i .Therefore the condition r Baumgardt et al. (2004a) also cited the further condition that the total mass of the cusp should not be larger than the mass of the IMBH, if the latter is heavier than a few percent of the cluster mass.From the point-mass model we obtain M NL  3500 M e at the 3σ level, i.e., about 1.4% of the cluster mass.Under the reasonable assumption that the true amount of nonluminous matter cannot be larger than twice the aforementioned upper limit, this constraint is not necessary in this case.A visual inspection of Figure 5 in Ferraro et al. (2003) led us to place an upper limit on r i .In fact, if the cusp is present, it should leave a signature in the optically observed radial profile of the cluster.In our inspection of the aforementioned figure, we could not identify any deviation from the King profile that describes the data at angular radii up to 10″.We thus deduce that, if indeed a cusp is present in the core, its radius of influence cannot be larger than the x-coordinate x 1 of the first data point in the plot, which we extrapolated to be = x log arcseconds 0.1 10 1 ( ) , i.e., 0.025 pc at the cluster distance.We conservatively assumed that the cusp can extend up to the x-coordinate x 2 of the second data point, which we extrapolated to be = x log arcseconds 0.6 10 2 ( ) , i.e., 4 0, which corresponds to 0.08 pc at the cluster distance.This value is roughly equal to the projected distance of PSRB from the clusterʼs center of gravity, the core pulsar with the second highest projected distance.One can easily verify that in the cases of PSRB, PSRD, and PSRE, the parameter space that results from all the considerations above also contains a subspace where these pulsars are located inside the cusp, thus ensuring that these three objects can be used as probes of its structure.
The resulting clusterʼs mass profile M(r), including the IMBH, is then: where the King profile that describes the mass distribution at r r i is the inner King profile from the studies by Ferraro et al. (2003;see Section 4.3).In this model the still unconstrained structural parameter is M BH , hence s = (M BH ), and we assumed its prior P M BH ˜( ) to be flat and nonzero in the conservative range 0 M BH 10 4 M e .We also included σ 0 among the parameters of our probability space, with a Gaussian prior whose mean and standard deviation are its measured value and 1σ uncertainty, respectively.We obtained space.The black line in the color plot marks the boundary between the two regions where the observational constraint on the IMBH radius of influence is (left side) and is not (right side) satisfied.We immediately see that the points with the highest probability density lie at the boundary of the allowed region, and not inside it, and this leads us to put into question the physical validity of our result, i.e., if there is indeed a ∼2.7 × 10 3 M e IMBH at the center of the cluster, or whether it is simply a mathematical solution for our problem with no real physical meaning.Figure 3 also displays the posterior probability of the M BHσ 0 space, but in the case where the constraint on the IMBH radius of influence is not applied.The points with the highest probability are located in the not allowed region, and those points, that show in Figure 2 the highest probability, have now a probability not negligibly lower than the maximum.Quantitatively speaking, the probability that the true solution of the problem lies in the allowed region, i.e., that the IMBH scenario has a real physical meaning, is ∼0.25 only.This probability is not low enough for an immediate ruling out of this model, nor high enough to consider physically certain the presence of an IMBH in the core of NGC 6752.We could in principle relax the upper limit on the radius of influence, allowing r i to be larger than the 4 0 limit assumed above.This would shift, in Figure 3, the diagonal line toward larger values for M BH , thus increasing the probability that the IMBH scenario has a real physical meaning.But a cusp with such a large radius would already have been detected in the analysis by Ferraro et al. (2003), at odds with their results.
If we, instead, put a more stringent upper limit on r i , assuming that its angular size at the cluster distance is not larger than the first data point of Figure 5 in Ferraro et al. (2003), i.e., r i x 1 = 1 25, further considerations arise against the physical validity of the IMBH scenario.It must be noted that, with such a tight constraint on r i , the cusp radius would be smaller than the projected distance of PSRE from the clusterʼs gravity center.We recall here that PSRE is the core pulsar with the smallest angular distance from the cluster center.This means that none of the core pulsars would be located inside the cusp, hence none of them could probe the cusp structure.As a consequence, the IMBH model would be totally indistinguishable from the point-mass model, and it should return exactly the same results for the overall amount of the nonluminous mass.It is easy to check whether this was the case, namely if we could obtain an overall amount for the nonluminous of at least 2.66 × 10 3 M e (3σ lower limit, see Section 4.4), by applying our methods with this new tighter constraint on r i .If we conservatively impose that the central velocity dispersion σ 0 lies in the range 6.5-10.5 km s −1 , i.e., within 10σ from its measured value of 8.5 km s −1 (see Table 3), the overall a priori predicted amount of the nonluminous mass can be at most 1.3 × 10 3 M e , in clear contradiction with the results given by the point-mass model.The IMBH model can, again a priori, predict M NL 2.66 × 10 3 M e , only if the clusterʼs central Posterior probability for the M BH -σ 0 space resulting from the analysis aimed to constrain the possible presence of an IMBH in the core of NGC 6752.The horizontal and vertical scales have been zoomed in the ranges 2 × 10 3 M e -4 × 10 3 M e and 8 km s −1 -9 km s −1 , respectively, for a better inspection of the relevant portion of each plot.The diagonal black line in the colored plot delimits the observational constraint on the size of the IMBH influence radius.The scale for the 2D color map, displayed in the vertical left panel, is normalized to the maximum value of the 2D probability distribution.Units for all probability distributions are arbitrary.
velocity dispersion is at least 15.1 km s −1 , which is too large with respect to the one inferred from the optical observations.
All of these considerations point against the presence in the core of the cluster of an IMBH, either with a mass equal to or higher than few thousands solar masses, or anyway massive enough to be responsible for the presence of the extra mass that explains the measured derivatives of the spin period for the core pulsars.The results from the point-mass model firmly state that an amount of nonluminous mass of at least ∼3 × 10 3 M e is necessary to justify the observed accelerations and jerks acting on the core pulsars, but it is very unlikely that the IMBH model is able to provide this amount of mass, consistently with all other optically observed features of NGC 6752.Most likely, the only possible scenario remaining invokes the presence in the cluster core of a population of much lighter nonluminous objects, whose nature and spatial distribution cannot be investigated with the present data, and which are so far undetected in optical observations and not yet predicted by models and simulations published by other authors.

The Central Mass-to-light Ratio of NGC 6752
The mass-to-light ratio (   ) in GCs provides a useful link between the directly measurable emitted light and the indirectly determinable distribution of mass, thus adding elements in the investigation of the structure of this kind of stellar associations, and probing models obtained by theoretical calculations and numerical simulations.
Following the method reported in Paper I, we revisited the mass-to-light ratio argument for PSRB, D and E, also including PSRF.Equation (1) in Paper I reads: where θ ⊥ is the projected angular offset of the pulsar with respect to the cluster center, q < âl,max ( ) is the maximum acceleration, due to the cluster potential that an object can experience at the angular offset θ ⊥ , M cyl (<θ ⊥ ) is the cluster mass enclosed in a cylinder of radius θ ⊥ whose main axis is parallel to the line of sight and passes through the cluster center, Σ V (<θ ⊥ ) is the clusterʼs average surface brightness within the angular radius θ ⊥ from its center, and L V,e is the luminosity of the Sun in the V-band.Using the value Σ V = 4.025 × 10 4 L e pc −2 , assumed constant according to the work reported in Noyola & Gebhardt (2006), we obtained >    4.6, 5.2, 4.6, and 4.3 (here and hereafter in this section in solar units for PSRB, PSRD, PSRE, and PSRF, in the given order).The discrepancies between these values and the ones reported in Paper I for the first three pulsars, namely     9 for PSRB and PSRE, and     13 for PSRD, are due to our use of more recent determinations of the surface  2, but also considering those points in the M BH -σ 0 space, whose corresponding value for the IMBH influence radius is larger than the upper limit we imposed on r i (see Section 4.5).Labels have been added for marking and easily identifying in which region our condition on r i is obeyed ("ALLOWED") or violated ("NOT ALLOWED").In both the uppermost and rightmost panels, the black solid and red dashed lines are the posterior probability distributions that are obtained by considering, in the marginalization, the allowed region only and the entire M BH -σ 0 space, respectively.All 1D posteriors have been plotted before being normalized after the marginalization, in order to highlight their relative height.
brightness profile and the coordinates for the cluster center (see Table 3).The position of the cluster center of gravity, assumed by the authors of Paper I, also led them to obtain the much larger lower limit of 1.3 × 10 4 M e for the mass responsible for the accelerations acting on PSRB, PSRD, and PSRE, than the one we obtained in our analysis (see Section 4).
We also derived the    values by considering the luminous mass distribution for the cluster reported and commented in Section 4. We integrated the double King density profile (see Section 4.3) along cylinders of radii equal to the angular separation of each core pulsar, and parallel to the line of sight.At first, we considered the luminous mass only, and we found =    2.0 for PSRB, PSRD, and PSRE, and 1.9 for PSRF.These values are substantially smaller than those we obtained by applying Equation (17).We then added a mass of 2680M e in the core, namely the 3σ lower limit on the nonluminous mass that results by applying the point-mass model, and we obtained    = 5. 65, 8.15, 9.31, and 4.35, thus confirming that Equation (17) may give a lower limit that is too conservative for    , hence for the amount of mass responsible for the acceleration acting on a pulsar in a GC.
It is worth mentioning a plot of the radial profile of    for NGC 6752, recently published in the online catalog of the Galactic GCs (see reference "d" in Table 3).This profile was obtained by Baumgardt (2017) by fitting the observed surface brightness profile against N-body simulations of the evolution of GCs, including the formation of nonluminous objects.We extrapolated from this plot the    values at the projected distance r ⊥ of each core pulsar, obtaining    = 3.1, 3.2, 3.3, and 2.6.These values are larger than the ones we obtained by considering the luminous mass only, yet still smaller than the ones implied by the  P P ratios and, consequently, than the ones we obtained considering our conservative lower limit on the amount of nonluminous matter.These discrepancies clearly are nonnegligible, and seem to mean that the methods used in Baumgardt (2017) still lead to an underestimation of the mass in the core of this cluster, in the form of objects that are not detectable in optical observations.

The Sources of Uncertainty in Our Results
In our Bayesian analysis we encountered three sources of uncertainty for our results.The first one is our poor knowledge of the intrinsic  P, for the core pulsars, which we addressed with a correction of the uncertainty in the measured accelerations.The second one resides in the treatment of the contribution due to the nearest neighbor to the time derivatives of the accelerations, a contribution that, as we already commented on, can be only statistically addressed.The third one is due to the uncertainty on the coordinates of the center of gravity of NGC 6752, which we have not taken into account so far.In this section we discuss the impact of these three sources of uncertainty on our results.

The Intrinsic Spindown
We treated our poor knowledge of the intrinsic spindown of the core pulsars in terms of an uncertainty to be added in quadrature to the one on all other terms on the right-hand side of Equation (3).As we already discussed in Section 4.2, we based our quantification of this additional term on the observed distribution of the measured  P for the MSPs in the Galactic field.Table 4 details all contributions to the uncertainty in the acceleration due to the GC, i.e., the uncertainty on the value one obtains by calculating a l from Equation (3), in terms of the 1σ fractional uncertainty.The term due to the intrinsic spindown dominates the uncertainty on a l by 2 orders of magnitude, with respect to the Galactic acceleration and the Shklovskii effect, and by 3 orders of magnitude with respect to the uncertainty on the measured  P. Nevertheless, the resulting total fractional uncertainty on a l is just a few percent, with an average of 3.9%.The resulting fractional uncertainty on the nonluminous mass can be estimated as follows.Equation (5) implies that the fractional uncertainty on a l is the same as on M TOT (r), being the latter the total amount of mass responsible for the observed acceleration, i.e., the sum of the amount of luminous mass M L with the amount of nonluminous mass M NL .A reasonable reference value for the amount of the luminous mass can be estimated by calculating the mass predicted by the King profile we used to describe the luminous mass in the core of the cluster (see Section 4.3), within a radius 〈r 3D 〉 equal to the average 3D distance of the core pulsars from the center of gravity.We estimated it as á ñ º á ñ = r r r 3 2 0.663 3D ,PSR c , where 〈r ⊥, PSR 〉 is the average angular displacement of the core pulsars from the coordinates of the clusterʼs gravity center.We thus obtained M L (〈r 3D 〉) = 460M e .By using the value M NL = 3050M e , namely the most probable value for the amount of nonluminous mass we obtained by applying the point-mass model, it results in M TOT (〈r 3D 〉) = 3510M e , whose 3.9% error is 137M e .Because we considered known the distribution of the luminous mass, this value represents the contribution  s M P , NL to the uncertainty on the nonluminous mass that is induced from the poor knowledge of the intrinsic spindown of the core pulsars., present the same quantities as the corresponding ones in Table 4, but for the case of the jerk instead of the acceleration, and P ̈instead of  P. The contribution to the uncertainty due to the intrinsic P ̈is not present since, as we already showed in Section 4, this quantity is several orders of magnitude lower than the uncertainty on P meas ̈, hence it can be completely neglected.The parameter  s a NN quantifies the uncertainty on the contribution to the jerk due to the nearest neighbor star.We defined it so that the probability to have

The Nearest Neighbor Contribution to the Jerks
| is equal to the probability that a generic Gaussian distributed quantity x is within ±1σ x from its mean x m , namely: ´-- ( ), where ρ K (r) is given by Equation (13).In this case the contribution due to the nearest neighbor dominates the overall uncertainty in  a l , while the uncertainty on , where η = 0.5404 is the average of the values in the fourth column of Table 5, the contribution on the uncertainty on M NL , due to the poor knowledge of  a NN only, results: It must be noted that the explicit expressions for  a l and  ¶ ¶ a M l NL also depend on the orientation, with respect to the observer, of both the pulsar position r and velocity v in the cluster's center reference frame.Because these dependencies cannot be neglected, we averaged Equation (20) over all possible orientations of the two vectors r and v, thus obtaining: We used again the results from the point-mass model (see Section 4.4), and the reference position = á ñ r r 3D , as above, thus obtaining must consider that the overall probability distribution for M NL is the product of the two distributions that result by separately considering the acceleration and jerk.One can easily see that the product of two peaked distributions, that are peaked at values that are close to each other but with different widths of the peak, results in another function that is peaked at the average value of the peaks of the two factors, with a characteristic width which is smaller than the one of both factors.Nevertheless, the value obtained for  s á ñ M a , NL NN seems so large that one may deduce that the probability distribution obtained by considering the jerks only might not be able to tighten the uncertainties on the nonluminous mass that is obtained by considering the accelerations only.For this reason we repeated the analysis reported in Section 4.4, but considering only the accelerations.We at the 2σ and 3σ level, respectively).The best value is consistent, within the uncertainties, with the one reported in Section 4.4, but in this case, the uncertainties are notably larger than the ones obtained by including the jerks.In practice, the amount of nonluminous mass is mainly determined by the accelerations, while the jerks allow a nonnegligible tightening of the uncertainties.

The Cluster's Gravity Center and the Final Estimate of the Amount of Nonluminous Mass
All quantitative results, reported and commented so far, rely on the assumed position for the cluster's gravity center, which we held fixed at the values indicated in Table 3 and hereafter referred to as the measured coordinates, whose 1σ uncertainties are 0 5 in both directions.Given their nonnegligible uncertainties, we investigated whether our estimation of the nonluminous mass might change if the center of gravity were not placed at the currently bona fide coordinates.We considered positions for the gravity center with an angular displacement r off 1″.For each position, we recomputed the amount of nonluminous matter by applying the point-mass model.In several cases, the resulting probability distribution showed more than one single peak for M NL .For this reason we consider in this discussion the mean mass M NL,mean (α, δ) so defined: where P(M NL , α, δ) is the probability density distribution obtained for the gravity center located at the celestial coordinates (α, δ).We obtained values in the range 2.87 × 10 3 M e  M NL,mean  7.54 × 10 3 M e , with mean 〈M NL,mean 〉 = 3.9 × 10 3 M e , standard deviation ˜. Figure 4 (panel (a)) displays a histogram of the obtained values, with a bin width of 200M e .There is a clear peak around 3000M e , in agreement with our previous result, but also a nonnegligible tail extending up to ∼7600M e .Figure 4 also contains information about the r off values, in steps of 0 25, from which one deduces that the tail above 5000M e receives contributions from positions at r off 0 5, and that the peak at 3000M e is mainly due to positions with an offset in the range 0 25 r off 0 75, with also for propagating the uncertainties between two quantities q and q 1 related to each other by q = f (q 1 , q 2 ,K,q N ), in the case where all other quantities q 2 ,K,q N are kept fixed.
relevant contributions from nearer and farther points.Such a distribution may be a consequence of our division of the exploited area in annuli, whose areas are in the ratio 1:3:5:7.In order to have a better insight into the dependence of the amount of nonluminous mass on the position of the cluster center of gravity, we plotted in Figure 5 a color map of the obtained values Figure 5. Amount of the nonluminous mass in the core of NGC 6752 as a function of the position of the cluster center of gravity.Panel (a): region containing the explored positions for the cluster center of gravity (small colored circle), compared to the position of the core pulsars (black bullets), and the core radius (blue dashed circle).Panel (b): a color map for the amount of the nonluminous mass in the core of the cluster as a function of the position of the cluster center of gravity.For each exploited position, the plotted value is the mean value of the posterior probability distribution (see Section 6.3).Dotted-dashed circles delimit the regions whose points have offsets r off 0.5σ, 0.5σ < r off 1.0σ, 1.0σ < r off 1.5σ, 1.5σ < r off 2.0σ from the innermost circle to the outermost annulus, respectively.In both panels, coordinates are centered at the measured position for the cluster center of gravity (RAJ = 19:10:52.04,DECJ = −59:59:04.64,Ferraro et al. 2003).The reference color scale is displayed in the top horizontal band.
for M NL,mean as a function of the position of the center of gravity.In panel (a) the obtained color map is placed in the context of the cluster core, while panel (b) shows an insight into the region under consideration.The position of the measured coordinates lies close to the northern boundary of a clear large area where M NL,mean (α, δ) is up to ∼3500 M e .Its displacement in decl.and elongation in R.A. can be seen as related to the fact that this area covers decl.around the mean of the core pulsars.The highest values are mostly required for positions toward the northwest, where one finds the farthest points from PSRD and PSRF, which are the pulsars with the highest  P | | (PSRD) and distance (PSRF) from the measured coordinates.
A final insight can be obtained after calculating the overall probability distribution P(M NL ), by marginalizing P(M NL , α GC , δ GC ) over the position of the cluster center.We performed this sum by taking into account the probability that the true center of gravity of the cluster is located in the considered position (weighted case): where P(α GC ) and P(δ GC ) are Gaussian distributions with mean value α GC and δ GC , respectively, and standard deviation σ = 0 5 in both cases, and also in the case of a flat distribution for the probability density of both coordinates (unweighted case): The top and middle panels of Figure 6 display the obtained unweighted and weighted distributions, jointly with the distribution for the nonluminous mass we obtained in Section 4, i.e., in the case where the cluster center of gravity is located at the measured coordinates.Both P W and P U peak at the same value of 2.97 × 10 3 M e , and have similar lower uncertainties at 1σ (−0.23 × 10 3 M e and −0.21 × 10 3 M e ), at 2σ (−0.33 × 10 3 M e and −0.29 × 10 3 M e , respectively) and 3σ (−0.41 × 10 3 M e and −0.39 × 10 3 M e , respectively).The higher mass tails are clearly not negligible up to 8 × 10 3 M e but, as shown by the histograms in Figure 4 and panel (b) in Figure 5, an Figure 6.Combined unweighted (red line, top panel) and weighted (green line, middle panel) probability distributions for the nonluminous mass in the core of NGC 6752, determined as in Section 6.3.The probability distribution obtained using the coordinates of Ferraro et al. (2003) for the cluster center is also presented for immediate comparison (blue line, bottom panel, label measured center).Units for the probability distributions are arbitrary.amount of nonluminous mass of at least 5 × 10 3 M e is required only if the true gravity center of the cluster has an offset r off 0 5 from the measured coordinates, i.e., 1σ away from the position measured by Ferraro et al. (2003) at the 1σ level, with a 3σ lower limit of 2.56 × 10 3 M e .This is the most constrained value for M NL we can obtain from our data after considering all the sources of uncertainty, yet its errors are much larger than the one introduced by our poor knowledge of the intrinsic first and second-order time derivative of the spin period of the core pulsars.A decrease by a factor of 10 on the uncertainties on the gravity center coordinates is still not enough for the uncertainty on M NL to be equal to the one induced by the the poor knowledge of  P. Adopting a very conservative approach, we do not conclude with a definite value for M NL , since from both P W and P U we deduce a 3σ upper limit on M NL at least of 10 4 M e .At the same time, we can place the solid lower limit of M NL 2.56 × 10 3 M e on the amount of nonluminous mass in the core of NGC 6752.
The position of the cluster's center of gravity is thus the dominant source of uncertainty on the estimation of the nonluminous mass in the core of NGC 6752.If its coordinates were known with a precision comparable to the one of the position of the pulsars in the core, the poor knowledge of both  P and P ̈for the core pulsars would become the main sources of uncertainty.But in such a situation, the resulting uncertainty on M NL would be of only a few percent in a case like NGC 6752, thus providing a good test bed for modeling and simulations.

Summary
We presented a timing analysis, spanning about 21 yr of observations, of the isolated pulsars PSR J1910-5959B, PSR J1910-5959C, PSR J1910-5959D, PSR J1910-5959E, and PSR J1910-5959F in the GC NGC 6752.The measured spin period derivative for the pulsars in the core is dominated by the gravitational pull of the mass in the central regions, and we estimated the amount of mass that is required to fully justify these measurements.We found that at least 2.56 × 10 3 M e , not yet taken into account in models and/or simulations, are present under the form of nonluminous matter, within a sphere of radius not larger than the distance of the closest pulsar to the cluster center of gravity.We also explored the scenario where an IMBH resides at the center of gravity of this cluster.We found that the presence of an IMBH with a mass of at least a few thousand solar masses can explain the observed accelerations and jerks experienced by the core pulsars, but is highly, although not completely, incompatible with the observational constraints deductible by the works of other authors, and that all these constraints hardly allow this scenario to provide all the necessary extra mass.We thus inferred that this extra amount of matter is most likely in the form of a population of much lighter sub-or nonluminous objects, whose exact nature and distribution cannot be investigated with the present data.Once this extra mass is taken into account, the inferred mass-to-light ratio in the central regions is higher than that jointly deduced from optical observations and from simulations of the structure and evolution of this GC.We have also shown that the main source of uncertainty in our results is the error in the coordinates of the cluster's center of gravity, and that our poor knowledge of the intrinsic spindown of the core pulsars has a smaller impact, unless the precision on the coordinates of the cluster center were comparable to that of the core pulsars, as one can obtain from their timing.

Figure 1 .
Figure 1.Timing residuals vs. MJD.ToAs obtained with Parkes and MeerKAT are plotted with blue triangles and red squares, respectively.If not visible, error bars are smaller than the marking symbol.Units and scales are the same for all pulsars for an immediate comparison of their timing rms and noise level.Pulsar names and rms values are indicated near the bottom right corner of the dedicated panel.
r δ ≡ δ PSR − δ GC , and r l , namely the pulsar distance from the center along the line of sight; the relation = ) and v δ ≡ dr δ /dt = (μ δ,PSR − μ δ,GC )D, where in turn D is the cluster distance, m d a cos ,PSR PSR and μ δ,PSR are the proper motion components of a pulsar, while m d a cos ,GC )).The adopted values for the Solar motion in the Galaxy are Θ e = 240.5 ± 4.1 km s −1 and R e = 8.275 ± 0.034 kpc (Gravity Collaboration et al. 2021).At the clusterʼs distance and position (l g = 336.4929deg, b g = −25.628deg in Galactic coordinates) we obtain a l,MW = (+2.7 ± 0.6) × 10 −11 m s −2 .The Shklovskii effect gives a contribution = a −2 , where μ is the proper motion amplitude, for a pulsar at the distance of NGC 6752.As already seen above,  P intr can be as high as several 10 −20 s s −1 , a value that would give a ∼5% contribution on the observed accelerations.Therefore, we calculated a l,m as  we took into account our poor knowledge of the intrinsic spindown by adding in quadrature the quantity   s ´á ñ + c P P P ( the 2σ and 3σ levels, respectively), while the overall amount of the nonluminous mass, M NL = M BH + M CUSP − M King (r i ), now results in = the 2σ and 3σ levels, respectively).Figure2displays the triangular plot for the posterior probability of the M BH -σ 0

Figure 2 .
Figure2.Posterior probability for the M BH -σ 0 space resulting from the analysis aimed to constrain the possible presence of an IMBH in the core of NGC 6752.The horizontal and vertical scales have been zoomed in the ranges 2 × 10 3 M e -4 × 10 3 M e and 8 km s −1 -9 km s −1 , respectively, for a better inspection of the relevant portion of each plot.The diagonal black line in the colored plot delimits the observational constraint on the size of the IMBH influence radius.The scale for the 2D color map, displayed in the vertical left panel, is normalized to the maximum value of the 2D probability distribution.Units for all probability distributions are arbitrary.

Figure 3 .
Figure3.Same as Figure2, but also considering those points in the M BH -σ 0 space, whose corresponding value for the IMBH influence radius is larger than the upper limit we imposed on r i (see Section 4.5).Labels have been added for marking and easily identifying in which region our condition on r i is obeyed ("ALLOWED") or violated ("NOT ALLOWED").In both the uppermost and rightmost panels, the black solid and red dashed lines are the posterior probability distributions that are obtained by considering, in the marginalization, the allowed region only and the entire M BH -σ 0 space, respectively.All 1D posteriors have been plotted before being normalized after the marginalization, in order to highlight their relative height.
results by using in Equation (9) the distance á ñ = r r 0.663 3D c from the cluster center, as in Section 6.1, σ V = σ V,PSR = 8.0 km s −1 , namely the 1D velocity dispersion of the pulsars for the local velocity dispersion at = á ñ r r 3D (see Section 3), and the approximation r á obtained M NL = (3.105± 0.137) × 10 3 M e at the 1σ level (

Figure 4 .
Figure4.Histograms of the amount of nonluminous matter in the core of NGC 6752 for the explored positions of the cluster center of gravity.Panel (a): cumulative histogram: for each bin different colors represent the contribution at distances r off , from the coordinates of the cluster center of gravity byFerraro et al. (2003), r off < 0 25 (red), 0 25 r off < 0 5 (dark blue), 0 5 r off < 0 75 (orange), 0 75 r off < 1 0 (light blue), respectively.Panels (b) to (e): separate histograms for each annulus defined as above.

Table 1
List of the MeerKAT Observations of NGC 6752 Used for This Work The symbols stand for the following: f c , central frequency; Δf, observing bandwidth; CD DM: dispersion measure used for coherent dedispersion; N chan , number of frequency channels; N pol , number of Stokes parameters; t samp , sampling time; N ant , number of antennas.

Table 2
Measured and Derived Parameters for the Pulsars in NGC 6752 Ferraro et al. (2003;ours, minutes, and seconds, and units of decl.aredegrees,arcminutes, and arcseconds.bReferenceepochforboth the spin period and the position.cTheoffset of the pulsars is calculated with respect to the position of the clusterʼs center of gravity reported byFerraro et al. (2003;see Table3).

Table 3
Positional, Kinematical, and Structural Parameters for NGC 6752 Published in Literature and Used in This Work

Table 5
P meas ̈plays a nearly negligible role.The resulting contribution l NN

Table 4
Contributions to the Uncertainties on the Pulsar Accelerations due to the GC Potential

Table 5
Contributions to the Uncertainties on the Pulsar Jerks due to the GC Potential 20 Equation (19) is the application of the known rule s = q . If one extracts from the probability distribution P W (M NL ) a measure for M NL with its (asymmetric) uncertainties, one obtains = M NL