Prospects for probing small-scale dark matter models with pulsars around Sagittarius A*

Future observations with next-generation large-area radio telescopes are expected to discover radio pulsars (PSRs) closely orbiting around Sagittarius~A* (Sgr~A*), the supermassive black hole (SMBH) dwelling at our Galactic Center (GC). Such a system can provide a unique laboratory for testing General Relativity (GR), as well as the astrophysics around the GC. In this paper, we provide a numerical timing model for PSR-SMBH systems based on the post-Newtonian (PN) equation of motion, and use it to explore the prospects of measuring the black hole (BH) properties with pulsar timing. We further consider the perturbation caused by the dark matter (DM) distribution around Sgr~A*, and the possibility of constraining DM models with PSR-SMBH systems. Assuming a 5-year observation of a normal pulsar in an eccentric ($e=0.8$) orbit with an orbital period $P_b = 0.5\,$yr, we find that -- with weekly recorded times of arrival (TOAs) and a timing precision of 1 ms -- the power-law index of DM density distribution near the GC can be constrained to about 20%. Such a measurement is comparable to those measurements at the Galactic length scale but can reveal small-scale properties of the DM.


I. INTRODUCTION
Black holes (BHs) are among the most fantastic objects in the Universe.They are holding important clues to some open questions in fundamental physics, concerning the curved spacetime, also possibly the connection of gravitation to the quantum world [1][2][3][4].From astrophysical observations, it is believed that supermassive BHs (SMBHs) exist at the center of most of massive galaxies [5].They provide us with precious opportunities to probe new physics beyond the current paradigm [6,7].The SMBH residing in our Galactic Center (GC), Sgr A*, has been confirmed with observations of S-star orbits as well as the image of its shadow [8][9][10][11].As a BH with mass around 4.3 × 10 6 M ⊙ and at a distance to the Solar System of about 8 kpc [9], Sgr A* has the largest mass to distance ratio among the known BHs.Thus it is an ideal laboratory for precision BH physics.
In the General Relativity (GR), it is known that for an isolated stationary BH, the spacetime around it is totally determined by three parameters: the BH mass (M • ), spin (S • ) and electric charge.It is the so-called no-hair theorem [12][13][14].All uncharged BHs in GR, which are usually considered in astrophysics, satisfy the Kerr solution.As a result of the nohair theorem, all the higher-order multiple moments of a Kerr BH can be expressed via its mass and spin [15].In particular, there is a relation between the BH's dimensionless spin parameter, χ • ≡ cS • /GM 2  • , and its dimensionless quadrupole moment, q • ≡ c 4 Q/G 2 M 3 • , via [16], Therefore an independent measurement of a BH's mass, * lshao@pku.edu.cnspin and quadrupole can give a direct test of GR (see e.g. Refs. [17][18][19][20][21][22][23]).
The mass of Sgr A* can be measured with monitoring gases or stars orbiting it in close orbits [8,9].Observations of the S2-star's orbit and its spectroscopy in the last three decades have not only provided a measurement of the BH mass but also given a clear evidence of the Schwarzschild precession [24], which sets an upper limit on the extended mass inside the S2's apocenter of about 3000 M ⊙ [25,26].However, due to the large orbital radius of S2 and the complex environment near the GC, the spin or quadrupole of Sgr A* may be hard to be measured with near-future S2 observations [27].Another effort that has been done is resolving the shadow of Sgr A* [10], which is a test of the BH metric at the scale of several Schwarzschild radii.Using the size and shape of the shadow image, the first set of results from the Event Horizon Telescope Collaboration gave a consistent constraint on the mass of Sgr A* with that from S-stars [10,28] and constrained BH alternatives [11,29,30].
A powerful tool of measuring the properties of the BH spacetime around Sgr A* is observing a pulsar orbiting it in a close orbit [18,20,31,32].Taking the advantage of the high accuracy of the pulsar timing technology, it is expected to detect the spin and quadrupole moment of Sgr A* with proper pulsars [18,21,31].Searches of pulsars in the GC have been carried out several times over the last few decades but no pulsar within the inner parsec has been found yet [33][34][35][36][37][38].Although observational evidence and theoretical model both indicate that there could be a number of neutron stars in the GC [35], the large dispersion measures and highly turbulent interstellar medium in the GC region make the detection nearly impossible at the typical low-frequencies [39].The observational sensitivity at higher frequencies are limited by the steep spectrum of pulsar emission.Nevertheless, future observations with the next-generation telescopes, such as the Square Kilometre Array (SKA) and the next-generation Very Large Array (ngVLA), are still hopeful of finding those pulsars [18,21].Such a discovery could open a new avenue of testing gravity [17,19].
Because of the potential for breakthrough, it is crucial to develop a pulsar timing model for PSR-SMBH systems.Several pioneering works [18,31] have been done based on an analytic solution that consistently includes the periodic spin effects derived by Wex [40], which can be regarded as an elegant extension of the widely used Damour-Deruelle (DD) timing model [41] for this particular situation.However, it is in general hard to analytically include both the spin and quadrupole effects simultaneously, and no elegant solution has be found yet.Nevertheless, the influences caused by the quadrupole effect can be included as perturbations to the pulsar coordinate position as well as for an additional secular precession [18].In this paper, we are going to develop a preliminary pulsar timing model based on the direct numerical integration of the post-Newtonian (PN) equation of motion that includes the spin-orbital coupling and quadrupole interaction.Comparing to applying a fully GR calculation based on the Kerr metric [20], the PN formalism allows us to treat the BH's spin and quadrupole moment as independent variables, which is favored when testing alternative gravity theories.The numerical method is also more flexible when considering new contributions from other perturbations.Different from the analytic approach, which describes the system with parameters represented for separate effects, such as the periastron advance parameter, k, and the deformation parameters of the orbit, δ r and δ θ , in the DD model [41], the numerical method uses the true physical parameters that directly related to the system, which is more convenient when further extending the model.
As we have mentioned before, the complex environment near the GC may spoil the PSR-SMBH system.The angular momentum and quadrupole moment contributed by the stellar cluster can obscure the signal caused by the BH's spin and quadrupole, which requires the orbit of the pulsar to be inside 0.1-1 mpc in order to meaningfully measure the BH's properties [27,31].Moreover, a high fraction of objects with mass about 10 M ⊙ may make the tests of gravity problematic at all radii [27].A conservative method is to only use the timing data that the pulsar is around the periastron, where the SMBH dominates the orbital dynamics of the pulsar [31].Another possible approach is to include the main contribution of the external perturbations as part of the timing model and do parameter estimation simultaneously.Similar idea has been used to constrain the existence of intermediate-mass BH in the GC with the S2 orbit [42].
In this paper we consider another kind of important perturbation in PSR-SMBH systems, that is the dark matter (DM) distribution around the SMBH.In the standard Lambda Cold Dark Matter (ΛCDM) model, the galaxies are formed inside DM halos, which consist of DM particles/fields that account for about 27% of the mass-energy at present universe but with a yet unknown nature [43].While the DM distributions at large scales can be measured with galactic rotation curves or gravitational lensing, the small scale structures have seldom been constrained by observations [44].Cold-DM-only simulation showed that the DM density has a cusp ∝ 1/r in the halo center, which is a known feature of the so-called Navarro-Frenk-White (NFW) profile [45].However, baryonic processes like the adiabatic growth of the central BH or supernova feedback can modify the profile, while some other possible DM models, such as the ultralight DM model, also predict different central density distributions [44].Observations of the galactic rotation curve have shown both evidences for cusp and core-like distributions [46].A cuspy density profile around the GC may contribute enough DM mass inside the pulsar obit that can be detected by pulsar timing observation, while satisfies the constraints from the S2 observation [26,[47][48][49].Recent work using the high-precision timing results of millisecond binary pulsars to directly measure the Galactic acceleration and derive fundamental Galactic parameters also shows the potential of detecting the DM distribution with pulsar timing technique [50].In this paper, we explore the prospects of constraining DM models via timing a pulsar around Sgr A* by extending our timing model to include the DM contribution.We consider the DM perturbation as a spherical mass distribution at the leading order and ignore the triaxial deformation of the DM that may exist in the CDM model [51].
The remaining part of this paper is organized as follows.In Sec.II, we describe the orbital dynamics of a PSR-SMBH system based on the PN equation of motion.Section III presents the basic concept of pulsar timing as well as the various effects we take into consideration in our timing model.We give an inverse timing formula of this numerical timing model in Sec.IV and use it in Sec.V to get the expected measurement precision of the BH properties.In Sec.VI, we extend our timing model to include the effects of DM distribution around the GC and show the parameter estimation results obtained from the extended timing model.Finally, we conclude in Sec.VII.

II. ORBITAL DYNAMICS
Differently from the Newtonian case, the two-body problem in GR has no general analytic solution.One can obtain the approximate equations of motion for well separated systems through the so-called PN expansion [52].We invoke PN equation of motion for the two-body orbital dynamics [53][54][55], where r ≡ r PSR − r • is the relative coordinate position vector in the harmonic gauge, and t is the coordinate time.By writing down the above equation, we have restricted ourselves to the case that only the BH is spinning, as the pulsar spin has, in general, negligible effects on the orbital motion.The Newtonian acceleration in Eq. (3) reads rN = −GM n/r 2 , with M ≡ M • + m PSR , r ≡ |r|, and n ≡ r/r.Besides the Newtonian term, other terms in Eq. (3) represent contributions from higher-order PN terms, spin-orbit coupling, quadrupolar effects to the orbit, and so on.For our purpose, we consider the following terms, where ṙ ≡ dr/dt, v ≡ dr/dt, v ≡ |v|, η ≡ m PSR M • /M 2 , and ŝ is the unit vector pointing along the BH spin.In particular, in our simulation the acceleration at 2PN order, r2PN , which is smaller by O v 2 /c 2 than r1PN , and the acceleration at 2.5PN order, which leads to the emission of gravitational waves, are not included [55].The acceleration caused by the BH's quadrupole is numerically at the 2PN order.For consistency we should include at least the r2PN term, but our treatment could be enough for the purpose of forecasting the precision of measuring the BH properties.In measuring the parameters from real data, more terms are in demand.We leave this for future studies.
In addition, because the mass ratio m PSR /M • < 10 −6 , we will ignore the mass of the pulsar for the moment.By doing this, we do not consider the back-reaction of the orbit to the spin of the BH, which is proportional to η [53]; thus the spin of the BH stays constant in our simulation and the pulsar is moving like a test particle in the BH's spacetime.Also, we will have a constant position for the BH that, ṙ• = 0 and v = v PSR .With η ≃ 0, Eq. ( 3) simplifies to, Our approximations are all well justified if our purpose is only to assert the precision of measuring the mass, spin, and quadrupole of the Sgr A* via timing a pulsar around it.When one wants to conduct a practical timing model to fit the real time-of-arrival (TOA) data from such a system, depending on the orbital characteristics, one needs to reconsider these approximations and to take into account all terms that would lead to timing residuals larger than the noises in observation.In general, by including higher-order terms one may have better power to break degeneracy among parameters due to the varieties that would be introduced by these terms.In this regard, our treatment is on the conservative side.The same arguments apply to the treatment of the Einstein delay and the Shapiro delay in Sec.III.We wish to track down the effects from higher-order terms in future investigation.
The orbit of a pulsar around Sgr A* is described by the orbital period, P b , the eccentricity, e, and various angles that determine the orientaion; see Fig. 1 for definition of these angles.A Keplerian description of the orbit is understood to be an approximation to the instantaneous motion of the pulsar.In our simulation, the orbital elements are updated according to the instantaneous position, r, and instantaneous velocity, v.For example, the eccentricity, rigorously speaking, should be a function of time, e(t), due to various non-Newtonian acceleration terms in Eq. ( 3).If not mentioned explicitly, the values of orbital elements, for example in Eq. ( 5), refer to the reference time t = 0 in the simulation.
As a fiducial case, we study a PSR-Sgr A* system with the following parameters, The angles are quasi-randomly chosen, to avoid special orientations that might render the parameter-estimation problem  5) for more details.
degenerate. 1 We integrate the orbit for t total = 5 yr, that contains t total /P b ≃ 10 orbits.TOAs are extracted weekly and uniformly in time.For 5 years, we simulate 260 TOAs in total.In the simulation, the longitude of the ascending node is fixed to Ω = 0.This choice is generic, because in this case the rotation of the pulsar orbit around K0 does not affect TOAs if we ignore the proper motion of the Sgr A* [57].
The orbits with Eq. ( 5) are illustrated in Fig. 2. It is seen in the figure the precession of the orbits due to the PN acceleration r1PN and the frame-dragging effects from the BH spin [58,59].

III. PULSAR TIMING
In pulsar timing, the TOAs of pulsar pulses at the telescopes are connected to the proper rotation numbers, N, of the pulsar [41,57,[60][61][62].The timing model incorporates the various effects in the orbital motion and radiation propagation, and thus from the observed TOAs one can extract the underlying physical parameters.For the pulsar's proper rotation, we assume, where ν ≡ 1/P is the pulsar's spin frequency, ν is the time derivative of ν, and T is the proper time of the pulsar [41].Higher-order time derivatives are easy to include when needed.At the lowest order, the proper time of the pulsar is connected to the coordinate time t via, 1 The value of λ was inspired by Psaltis et al. [56].
Integrating the above equation gives an Einstein delay that accounts for the gravitational redshift and special-relativistic time-dilation effects [41,60], To drop a term that is linear in time, one can redefine the pulsar spin [41] and Eq. ( 7) takes a form as Roughly speaking, the ⟨ •⟩ term that appears in the denominator means to average over the pulsar orbital period, so that ⟨dT/dt⟩ = 1 and there is no linear-in-time dependence term in the Einstein delay.Due to the spin and quadrupole effects, the denominator is no longer a constant.Here we use the 1PN approximation of the denominator, which accounts for almost all the linear dependence.For a Keplerian orbit, ∆ E is found to be [41,60], where u is the eccentric anomaly.In our simulation, in order to account for all the factors coming from the variation of the orbit, we integrate Eq. ( 9) to obtain the Einstein delay.
The orbital motion of the pulsar gives the geometric delay, called the Römer delay, which is simply, We also include the lowest-order propagation time delay caused by the curvature of the Sgr A*, called the 1PN Shapiro delay [60,63], Finally, we assume that TOAs are collected at an infinite distance to the Sgr A*.In doing so, we are ignoring various terms related to the proper motion of the Sgr A* and the motion of the Earth around the Sun, etc..After dropping the constant (infinite) term, one finally has [41], Notice that we are only considering the lowest-order terms in the Einstein delay and the Shapiro delay.In principle, higher-order terms can also be added; for examples, see Kopeikin [64] and Wex and Kopeikin [6] for the higher-order terms in the propagation delay of pulses.We suspect that the higher-order terms will further break the degeneracies in the parameter-estimation problem which, as we mentioned, renders our treatment conservative.We defer the investigation for a future study.
We give an illustration in Fig. 3 for these three time delays for the fiducial orbit defined in Eq. ( 5).These delays are extremely large compared with those of the binary pulsar systems we are currently regularly timing [18,65].Notice that for our fiducial orbit, the inclination is only π/5 = 36 • .Even with such a small inclination, the Shapiro delay is already numerous; this is also true even for face-on orbits [18].FIG. 3. Three time delays in pulsar timing, defined in Eq. ( 11), Eq. ( 8), and Eq. ( 12), for the fiducial orbit (5).The linear trend in time in ∆ E is removed; it can be absorbed into a redefinition of pulsar's spin.

IV. THE INVERSE TIMING FORMULA
In the standard procedure of parameter estimation in pulsar timing, for calculating the residuals one needs an inverse timing formula that calculates N from given TOAs and system's parameters Θ instead of the pulsar timing model described in Sec.III that calculates t TOA (N; Θ) [41].In principle, giving the system's parameters Θ and the observed TOA, t TOA , we can first integrate the pulsar's orbital motion and get a series of time delays that are related to the coordinate time: ∆ R (t i ) , ∆ S (t i ) , ∆ E (t i ).Then we can do interpolation for these series and solve the implicit equation, t TOA = t + ∆ R (t) + ∆ S (t), to get the coordinate pulse emission time, t.Finally we use T = t − ∆ E (t) and Eq. ( 6) to get the related N.But a problem rised in doing this procedure numerically.To have an acceptable precision in solving the implicit equation, the interpolation step needs a dense coverage, which makes the orbital integration slow.So here we propose a fast method to get the inverse timing formula.
Briefly speaking, we want to change the variable t in the differential equations ( 3) and (7) to t TOA that is directly related to the observation.
From Eq. ( 8) and Eq. ( 13), we have Taking a derivative on both sides and insert the explicit forms shown in Eq. ( 11) and Eq. ( 12), we get Combine this with Eq. (3) and Eq. ( 7), we can get the complete differential equations used in the inverse timing formula, which are where dt/dt TOA is the inverse of Eq. (15).Note that we also need to transform the initial conditions.From the conditions that at t = 0, the system has parameters Θ and ∆ E = 0, we can simply calculate t TOA 0 = ∆ R t=0 + ∆ S t=0 .The initial conditions now are at t TOA = t TOA 0 , and the system has parameters Θ and ∆ E = 0.
Integrating the above equations from t TOA 0 to t TOA , we can get the related pulsar proper time T by and solve N from Eq. ( 6) easily.We will use this method as the inverse timing formula in the following numerical calculation.

V. PARAMETER ESTIMATION
After all parameters, which we denote as Θ, are given, one obtains the pulsar rotation number N as a function of t TOA without any ambigiuity.The parameters include, where Θ • includes parameters of the BH in Eqs.(5a-5b), Θ orbit includes parameters of the orbit in Eqs.(5c-5d), Θ PSR includes parameters of the pulsar spin in Eq. ( 6), and Θ pert includes parameters from the perturbation.Therefore, we denote the pulsar rotation number as N Θ; t TOA .Θ pert is to be introduced in Sec.VI, and we omit it in this section.
Here in Θ • , we treat χ • and q • as independent variables.In GR, the no-hair theorem poses q • = −χ 2 • [3,12].Our treatment can be viewed as an expansion of the multipoles of the spacetime [16] with the three lowest-order moments, namely the monopole M • , the dipole S = χ • GM 2 /c, and the quadrupole Q = q • G 2 M 3 /c 4 .The post-Newtonian expansion in Eq. (3) provides the possibility to treat χ • and q • independently, thus providing the possibility to assess the precision in testing the no-hair theorem.
Assuming a Gaussian timing noise realization in observation, the probability that the true values of parameters being Θ is, where , and the summation is over the number of TOAs.In real data, the signal is contaminated with noises, and the simulation can take them into account by adding, say, Gaussian noises, as was done in Refs.[66,67].Nevertheless, it is statistically equivalent to use the noiseless templates, N Θ; t TOA in Eq. ( 22).The addition of random We assume that the pulsar spin frequency ν = 2 Hz and the timing precision σ TOA = 10 −3 s.Simulations include weekly observations over a 5-yr interval.
noises only shifts ln P Θ t TOA by a constant value statistically for the expecting value of an ensemble of noise realizations.The use of noiseless templates avoids the randomness in the realization of noises, and it was also adopted in other scientific studies; see e.g.Gaebel and Veitch [68] for simulations of parameter estimation with gravitational-wave waveform templates.As a standard procedure for parameter estimation [41,69], we estimate the measurement uncertainties of these parameters via the covariance matrix where L = − ln P Θ|t TOA is the log-likelihood function.We assume the timing precision σ TOA = 1 ms in our simulation, which is a relatively conservative estimate for future observations [18].Figure 4 shows the simulation results for the fiducial system configuration except for the pulsar orbital period which varies from 0.1 yr to 1.0 yr.The range of P b chosen here was inspired by Liu et al. [18].Previous studies suggest that ∼ 10 3 pulsars can be expected to be orbiting around Sgr A* with P b < 100 yr [35], while the innermost some of them could be in orbits as tight as ∼ 100 − 500 au from Sgr A* [70], which correspond to the orbital periods considered here.From Fig. 4, we can see that for orbits with P b ≲ 0.5 yr, the measurement of the spin and quadrupole parameters can be better than 1%, which is consistent with the results in Refs.[18,21].
We also investigate the influence of the orbital eccentricity and the results are presented in Fig. 5.It shows that the fractional precisions of the mass, spin, quadrupole parameters are nearly linear functions of the factor 1 − e 2 in the logarithmic scale.This can be explained intuitively by that the secular effects of mass, spin, and quadrupole are proportional to Fractional Precision Fractional precision for the BH mass, spin and quadrupole parameters as functions of the orbital eccentricity.The system parameters (except e) are the same as before while the orbital eccentricity changes from 0.1 to 0.9.1−e 2 −1 , 1−e 2 −3/2 , and 1−e 2 −2 respectively [6,53].From this figure we can also conclude that for obits with very high eccentricities, say e ≳ 0.8, the precision of the quadrupole determination can reach the precision of the spin determination or even better, which comes from the ∝ r −4 nature of the quadrupole interaction.
As shown in Refs.[18,20], there are leading-order degeneracies among spin parameters.The observable secular effects caused by the BH spin including the advance of the periastron ω and the change of the inclination angle i [6].However, even if one can separate the periastron advance caused by the Schwarzschild precession, these two secular effects still cannot fully determine the three spin parameters.Higher-order secular effects, for example, the change of the periastron advance rate ω, or periodic effects are needed when measuring the BH spin and its orientation [18].In Fig. 6 we plot the correlation between the spin parameters, where blue lines give the predicted leading order degeneracies determined by [20] where S = S ŝ is the spin vector of the BH and Â is the unit vector pointing to the periastron from Sgr A*.Here we have used the initial orbital elements in Eqs.(5a-5d) to calculate the leading-order degeneracies.The small discrepancies between the predicted directions and the major-axis of the contours may come from the change of the orbital elements and the higher-order contributions in the precessions or from the periodic effects.

VI. DARK MATTER PERTURBATION
In previous sections we have assumed the system to be sufficiently clean and estimated the potential of testing gravity using such PSR-SMBH systems.In reality, there may be various effects that could complicate or even spoil these tests.Such as the gravity perturbation from massive objects and surrounding mass distribution [27].One way to avoid the external perturbation is to focus on the timing features near the periastron, where the interaction with central BH is dominating [31].In this work we choose to extend our timing model to include the parameters that describe the perturbations and study the actual influences of specific kinds of perturbations.Here we consider the perturbations from a spherically distributed DM mini-halo around the GC.The DM profiles are introduced in Sec.VI A and Sec.VI B. In Sec.VI C we present the extended timing model and the parameter estimation results.

A. The Generalized NFW Profile
The standard ΛCDM model of the universe has been remarkably successful in explaining the evolution of the universe and the development of large-scale structures.The cold-DM-only simulation gives a nearly mass-independent DM halo density distribution, the so-called NFW profile [45] where R s is the scale radius and ρ 0 is a characteristic density.This profile has been widely used in describing galaxies' DM halos and fitting for the mass distribution in lensing observations.For our Milky Way, it gives R s ≃ 20 kpc, and ρ 0 r≈8 kpc = 0.4 GeV cm −3 at the location of our Solar System [71].
A feature of the NFW profile is the singular density profile at the center where ρ(r) ∝ r −1 , which is called a density cusp.However, observations of the rotation curve of lowsurface-brightness galaxies suggest that some galaxies have a finite density core which can form in some other DM models such as the self-interacting DM (SIDM) model [46,72].Although the cold DM model has been successfully examined by observations at large scales, the small scale properties of DM are still lack of constraints [44].Thus a measurement of the central profile of DM density may help us understand it better and constrain different DM models.Various baryonic processes such as adiabatic contraction or density fluctuations due to supernova feedback could have modified the DM profile [73,74], while recent simulations also suggest that SIDM models do not produce large differences in the inner structure of Milky-Way-mass galaxies in the presence of baryonic feedback effects [72].Nevertheless, to account for different model predictions, here we consider a generalized model, the gNFW model [75] The central behaviors of different models are characterized by the power-law index γ.Gondolo and Silk [76] pointed out that the adiabatic growth of the SMBH in the DM halo center will strongly modify the DM profile inside the radius of its gravitational influence, R h ≡ GM • /v 2 0 , which is about 1.7 pc for Sgr A* [77].For collisionless DM with a polytrope phase-space distribution, a DM spike with ρ sp (r) ∝ r −γ sp will form, with γ sp = (9 − 2γ)/(4 − γ).For a pulsar with orbital period P b ∼ 1 yr has a semi-major axis ∼ 1 mpc, which is at the region that is dominated by the central BH.It motivates us to take this effect into account.
Some DM models allow the DM particles to annihilate with themselves [76].If so, due to the very high DM densities in the halo center, it may act as a strong gamma-ray source [76], and the spike induced by the SMBH will enhance it further.On the other hand, the annihilation of DM particles will limit the maximum density of the DM spike, and produce a weak cusp in the center with ρ in (r) ∝ r −γ in , where γ in ≃ 0.5 [78].Such a weak annihilation cusp will form around a radius R in where the DM density reaches a critical density ρ ann .
Combining the above discussion, we consider the DM density profile as below [79], If there is no DM annihilation, the profile is where R sp is chosen to equal to the gravitational radius R h .Densities ρ sp (r) and ρ in (r) are determined by the continuous condition at R sp .For the model with annihilation, we assume ρ sp (R in ) = ρ in (R in ) = 2ρ ann to interpolate these two profiles.

B. The Einasto Profile
The Einasto profile is also a commonly used DM distribution in which the logarithmic density slope shows a powerlaw behavior [80][81][82], which cannot be characterized by the gNFW model we discussed before.It is argued that this model provides a better fit to the high-resolution N-body DM simulation and it gives a core-like central behavior [81,83].The DM density profile of this model can be written as [48] ρ where ρ 0 is the DM density at the scale radius r s and α is the inverse of the Einasto index.Fitting the rotation curve and globular cluster kinematics from Gaia data gives [84] ρ The analytic form of the DM spike induced by the adiabatic growth of the SMBH in the Einasto model is not available in the literature.However, one may estimate it with circularorbit approximation [48].Following Shen et al. [48], we estimate the DM spike via where M DM (r) is the mass of DM enclosed in the radius r and M tot (r) = M SMBH + M DM (r).Exact distribution of the DM spike in the Einasto model can be numerically calculated as shown in Ref. [48].Considering that we are only interested in the central behavior of the DM spike at the radius scale of mpc, which is much smaller than the scale radius r s , it is a good approximation to have ρ E (r) = ρ 0 e 2/α .Combining with above equations, one derives that the DM spike density in the Einasto model is ρ E, sp (r) ∝ r 9/4 , which is effectively to have γ = 0 in the gNFW model.It is also demonstrated that even the initial Einasto profile is much smaller than the NFW profile at the central region, the DM spike densities in these two model in fact can have very similar orders of magnitude [48].
Even though the Einasto profile itself cannot be desicribed by the gNFW profile, we conclude that the DM spike distribution in the Einasto model can be well characterized by the first expression of Eq. ( 29), and two models will give basically the same results.
In Fig. 7 we show the DM density profile and the DM mass enclosed by the radius r.We take the NFW model as an example and show the related models with spike and DM annihilation.Using the method of osculating elements [85], one can calculate the secular effects caused by the DM distribution.Due to the spherical symmetry, the only secular effect caused by the DM is the modified advance of periastron.The precession time scales of various effects in PSR-Sgr A* systems are shown in Fig. 8. From Fig. 7 and Fig. 8, we can see that the DM mass contributed by the initial NFW model is far below the precision of BH mass determination obtained before, which is about 10-100 M ⊙ , while for the DM model with annihilation, timing a pulsar with an orbital period P b ≲ 10 yr can only give the information of the weak cusp caused by the DM annihilation, described by γ in , which is not strongly related to the original DM profile.For a pulsar with a larger orbital period, in principle one may have a measurement of γ sp , but the complex environments and external perturbations will complicate this system.For DM model without annihilation, the spike structure induced by the BH can largely increase the DM density, thus timing a pulsar with a reasonable orbital period will give a constraint on γ sp , which is related to the underlying DM profile via γ sp = (9 − 2γ)/(4 − γ).

C. Prameter Estimation
Compared to the mass of the central BH, the DM mass is still very small (about 10 −5 M • for the spike model inside 1 mpc).Thus we may treat the DM as a perturbation, which means that we only consider the Newtonian gravity caused by the DM distribution and find the effects on the pulsar's orbit.To verify this assumption, as an example, we calculate the leading-order Shapiro time delay caused by the extended DM mass distribution in the spike model, as shown in Fig. 9, which is for a system with parameters in Eqs.(5a-5d).One can simply estimate from Fig. 3 that the Shapiro time delay caused by the DM should be at the order of 1 ms, which is consistent with the calculation, while the extended mass distribution slightly weakens the amplitude and broadens its shape.The Shapiro delay caused by the DM is still almost degenerate with the Shapiro delay of the BH and its value is smaller than the assumed timing precision.Even with a larger orbit that includes more DM mass, we can still ignore this contribution in our setting.
In Fig. 10 and Fig. 11 we present the timing residuals caused by the DM perturbation.We simulated a set of TOAs with the full timing model that includes the DM effects and fitted for the TOAs with the model described in Sec.III, which does not include the DM distribution.This is the case for real observations that in the beginning people usually use a simple timing model that accounts for those largest effects as a first step.Figure 10 shows the residuals before fitting, i.e. we use the true parameters Θ that are used in simulation but do not include the DM contributions.Although the true parameters in principle are unknow, this figure shows the cumulated timing residuals caused by the DM perturbation.There is a secular part in the residuals due to the periastron advance caused by the DM perturbation.Figure 11 shows the residuals after fitting.The secular part of the residuals is mainly absorbed by the BH spin parameter and only a quasi-periodic part is left.We plot the residuals for t from the second year to the third year of the total observation time span, which corresponds to two orbital periods of the pulsar.The residuals after fitting remain an amplitude of about 1 ms, which is close to the assumed timing precision, showing a possibility of estimating the DM parameters with pulsar timing.The timing residuals for the initial NFW model and the spike model with DM annihilation have similar behaviors as show in Fig. 10 and Fig. 11, but with a much smaller amplitude as discussed above.The observation of S2 orbit has set a limit on the extended mass inside the S2's apocenter, which is about 3000 M ⊙ [25,26], and the DM model we considered here gives a value that is consistent with this constraint.The extended mass inside the S2 orbit contributed by the DM in the spike model is ≲ 1000 M ⊙ , although the star clusters will also contribute to the extended mass in this scale [26].
As discussed before, we extend our timing model by adding a DM distribution related to the DM spike structure, which is a density profile with two parameters where ρ 0 now represents the central density of the DM spike and γ sp is related to the original DM profile, where we use the NFW profile with γ = 1 in our simulation.We fit the TOAs with the full model where the effects from DM are included.In Fig. 12 we show the results of parameter estimation.The input parameters of the BH and the pulsar are shown in Eqs.(5a-5d), and the DM model is the NFW model with spike.We still assume weekly observations over a 5-yr interval.We can see that timing a pulsar with an orbital period P b ≳ 0.5 yr and an orbital eccentricity e ∼ 0.8 can give a fractional measurement uncertainty of ∼ 1% in γ sp for the spike model, which is related to a ∼ 20% fractional precision in γ.This result is comparable to the constraints from fitting kinematic data of maser and other observations, which give γ = 0.79 ± 0.32 for our Galaxy [71].However, by timing a pulsar around Sgr A*, we can constrain the DM structure in the length scale of ∼ mpc, which is 6 orders of magnitude smaller than the length scale of other Galactic observations of such an investigation, which are typically done at ∼ kpc scales.
We also investigate the effects of the orbital eccentricity and the result is shown in Fig. 13.We can see that the measurement precisions of the DM parameters decrease fast as the orbital eccentricity becomes small.This is because for small orbital eccentricity, the effect caused by the DM distribution is strongly degenerate with the BH mass, providing only a mass monopole in the limit of e → 0. Only when the pulsar moves at different radii it can sense the DM's radial distribution.
In the above, we have used unrelated χ • and q • , although in GR they can be related via Eq.( 1).Another consideration is to detect the DM profile under the assumption that GR is correct, which corresponding to setting q • = −χ 2 • in the timing model.The constraints on the DM parameters obtained under this assumption are similar to what are shown in Fig. 12 and Fig. 13.This can be partly explained by Fig. 8.As the only secular effect caused by the DM distribution is the periastron advance, the leading-order degeneracy among the DM parameters and other parameters are dominated by this secular effect, which is similar to the discussions on the leading-order degeneracy among the spin parameters.Thus the quadrupole effect can only have a very small contribution compared to the others as shown in Fig. 8. Also it has been discussed by Heißel et al.
[26] that, the orbital features caused by the extended mass distribution are significant in the orbital section θ ∈ (5π/6, 7π/6), which is different from the spin or quadrupole effects that are most significant near the periastron.So setting q • = −χ 2 • only has small effect on the DM measurement.However, if the central object is totally different from GR, such as if it is a supermassive boson star that can have larger quadrupole moment, the conclusion may be different.

VII. CONCLUSIONS
In this work, we explore the prospects of constraining the BH properties and DM models from timing a pulsar around Sgr A*.We construct a timing model based on the numerical integration of the PN equation of motion, with leadingorder effects of various time delays being taken into account.Our simulations show that for a pulsar with an orbital period P b ≲ 0.5 yr and an orbital eccentricity e ∼ 0.8, a 5-yr observation with weekly recorded TOAs and a timing precision 1 ms, the BH spin and quadrupole parameters can be measured with a precision of 10 −2 or better, which is consistent with previous studies based on a semi-analytic approach [18,31] or a fully general-relativistic treatment [20].We need to emphasize that the timing model in this work only considers the leading-order effects, which is sufficient for estimating the measurement precision of the system parameters but not enough for applying to the real observations.Higher-order effects should be carefully considered and added into the timing model if the observational precision is high enough.We hope to further develop the timing model along this line in future studies.As a concrete application-which is hard to achieve with the semi-analytic approach or the fully general-relativistic treatment-by extending our timing model with a spherical DM perturbation, we investigate the measurement precision of the DM distribution at small scales.For DM models with a spike structure induced by the adiabatic growth of the central BH, timing a pulsar with an orbital period P b ≳ 0.5 yr, and an orbital eccentricity e ∼ 0.8 can provide a 1% measurement of the power-law index γ sp of the DM spike, which relates to a measurement of the underlying DM model with a 20% precision in the index γ.Such a precision is comparable with, but complementary to, the observations at the Galactic scale, usually done at ∼ kpc scales, which can be an equilibrium based kinematic analysis [71], or some extreme-precision time-series measurements of Galactic accelerations [86,87].As well known, measuring the DM distribution is in general easier in larger scales, where the DM can contribute a larger total mass.However, the various environmental effects can complicate the test in this situation.With future discovery of proper pulsars near the GC, pulsar timing observation with high precision will provide us a unique opportunity to explore the small-scale properties of the DM, and eventually lead to a more complete understanding of the origin of the DM.

FIG. 1 .
FIG. 1. Coordinate system and notations that are used for a pulsar orbit around Sgr A*.The frame Î0 , Ĵ0 , K0 has K0 pointing from the Earth towards Sgr A*, and Î0 , Ĵ0 constitutes the sky plane with Î0 pointing towards east and Ĵ0 towards north.The orientation of the pulsar orbit is determined by the inclination, i, the longitude of periastron, ω, and the longitude of ascending node, Ω.The Sgr A* sits at the origin, and its spin points towards ŝ ≡ (sin λ • cos η • , sin λ • sin η • , cos λ • ) in the coordinate system; in the figure λ • and η • are not shown for clarity.The position of the pulsar at the reference time t = 0 is r 0 , and it is described by the true anomaly θ 0 .

FIG. 4 .
FIG.4.Fractional precision for the BH mass, spin and quadrupole parameters as functions of the pulsar orbital period.The BH parameters and orbital parameters (except P b ) are shown in Eqs.(5a-5d).We assume that the pulsar spin frequency ν = 2 Hz and the timing precision σ TOA = 10 −3 s.Simulations include weekly observations over a 5-yr interval.

FIG. 6 .
FIG.6.Correlations between spin parameters for the fiducial system given in Eqs.(5a-5d).The contours represent the 68% and 95% confidence regions.True values are marked with blue dots and the blue lines show the theoretical leading-order degeneracies.

FIG. 7 .
FIG. 7. The DM density profile (upper) and the DM mass inside r (lower) for the NFW model, the NFW model with DM spike and the model with DM annihilation.

FIG. 8 .
FIG.8.Precession time scales as functions of the orbital period for various relativistic effects and the DM perturbation.

FIG. 9 .
FIG. 9.The Shapiro time delay caused by the DM distribution in the spike DM model which is based on the NFW profile.The system parameters are given in Eqs.(5a-5d).

FIG. 10 .
FIG.10.Timing residuals caused by the DM perturbation before fitting.We use the DM model with spike but no annihilation, and we set R sp = R h .The system parameters are given in Eqs.(5a-5d).

FIG. 11 .
FIG. 11.Same as Fig. 10, but after fitting for a model without DM.Most of the residuals in Fig. 10 are absorbed in other parameters.

FIG. 12 .
FIG. 12. Fractional precision for the BH mass, spin, and quadrupole parameters, as well as the DM parameters.Parameters for Sgr A* and the pulsar orbit are shown in Eqs.(5a-5d).The parameter ρ 0 and γ sp are the parameters for the NFW model with spike.

FIG. 13 .
FIG.13.Similar to Fig.12, but as a function of the orbital eccentricity.