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

Zexin Hu  Department of Astronomy, School of Physics, Peking University,
Beijing 100871, China  Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing 100871, China  Lijing Shao<EMAIL_ADDRESS>Kavli
Institute for Astronomy and Astrophysics, Peking University, Beijing 100871,
China  National Astronomical Observatories, Chinese Academy of Sciences,
Beijing 100012, China  Max-Planck-Institut für Radioastronomie, Auf dem Hügel
69, D-53121 Bonn, Germany  Fupeng Zhang  School of Physics and Materials
Science, Guangzhou University, Guangzhou 510006, China  Key Laboratory for
Astronomical Observation and Technology of Guangzhou, 510006 Guangzhou, China
Astronomy Science and Technology Research Laboratory of Department of
Education of Guangdong Province, Guangzhou 510006, China

###### Abstract

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 Hawking and Ellis (2023); Chandrasekhar (1983); Chrusciel
_et al._ (2012); Unruh and Wald (2017). From astrophysical observations, it is
believed that supermassive BHs (SMBHs) exist at the center of most of massive
galaxies McConnell and Ma (2013). They provide us with precious opportunities
to probe new physics beyond the current paradigm Wex and Kopeikin (1999); Will
(2008). 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 Ghez  _et al._ (2008); Genzel  _et al._ (2010); Akiyama  _et al._
(2022a, b). As a BH with mass around $4.3\times 10^{6}\,M_{\odot}$ and at a
distance to the Solar System of about 8 kpc Genzel  _et al._ (2010), 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_{\bullet}$), spin ($S_{\bullet}$) and electric charge. It is the so-
called no-hair theorem Israel (1967); Carter (1971); Robinson (1975). All
uncharged BHs in GR, which are usually considered in astrophysics, satisfy the
Kerr solution. As a result of the no-hair theorem, all the higher-order
multiple moments of a Kerr BH can be expressed via its mass and spin Hansen
(1974). In particular, there is a relation between the BH’s dimensionless spin
parameter, $\chi_{\bullet}\equiv cS_{\bullet}/GM_{\bullet}^{2}$, and its
dimensionless quadrupole moment, $q_{\bullet}\equiv
c^{4}Q/G^{2}M_{\bullet}^{3}$, via Thorne (1980),

$q_{\bullet}=-\chi^{2}_{\bullet}\,.$ (1)

Therefore an independent measurement of a BH’s mass, spin and quadrupole can
give a direct test of GR (see e.g. Refs. Kramer  _et al._ (2004); Liu  _et
al._ (2012); Shao  _et al._ (2015); Zhang and Saha (2017); Bower  _et al._
(2018); Weltman  _et al._ (2020); Eatough  _et al._ (2023)).

The mass of Sgr A* can be measured with monitoring gases or stars orbiting it
in close orbits Ghez  _et al._ (2008); Genzel  _et al._ (2010). 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 Abuter  _et al._ (2020), which sets an upper
limit on the extended mass inside the S2’s apocenter of about
$3000\,M_{\odot}$ Abuter  _et al._ (2022); Heißel  _et al._ (2022). 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 Merritt  _et al._ (2010). Another effort that has been done is
resolving the shadow of Sgr A* Akiyama  _et al._ (2022a), 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 Akiyama  _et al._ (2022a, c) and constrained BH
alternatives Akiyama  _et al._ (2022b); Xu  _et al._ (2023a, b).

A powerful tool of measuring the properties of the BH spacetime around Sgr A*
is observing a pulsar orbiting it in a close orbit Zhang and Saha (2017); Liu
_et al._ (2012); Psaltis  _et al._ (2016); Dong  _et al._ (2022). 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 Liu
_et al._ (2012); Psaltis  _et al._ (2016); Bower  _et al._ (2018). 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 Kramer  _et
al._ (2000); Eatough  _et al._ (2013); Wharton  _et al._ (2012); Siemion  _et
al._ (2013); Liu  _et al._ (2021); Torne  _et al._ (2023). Although
observational evidence and theoretical model both indicate that there could be
a number of neutron stars in the GC Wharton  _et al._ (2012), the large
dispersion measures and highly turbulent interstellar medium in the GC region
make the detection nearly impossible at the typical low-frequencies Cordes and
Lazio (2002). 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 Liu  _et al._ (2012); Bower  _et al._ (2018). Such a discovery
could open a new avenue of testing gravity Kramer  _et al._ (2004); Shao  _et
al._ (2015).

Because of the potential for breakthrough, it is crucial to develop a pulsar
timing model for PSR-SMBH systems. Several pioneering works Liu  _et al._
(2012); Psaltis  _et al._ (2016) have been done based on an analytic solution
that consistently includes the periodic spin effects derived by Wex (1995),
which can be regarded as an elegant extension of the widely used Damour-
Deruelle (DD) timing model Damour and Deruelle (1986) 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 Liu  _et al._ (2012). 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 Zhang and Saha
(2017), 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, $\delta_{r}$ and $\delta_{\theta}$, in the DD model
Damour and Deruelle (1986), 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 Psaltis  _et al._ (2016); Merritt
_et al._ (2010). Moreover, a high fraction of objects with mass about
$10\,M_{\odot}$ may make the tests of gravity problematic at all radii Merritt
_et al._ (2010). 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 Psaltis  _et al._ (2016). 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 Straub  _et al._ (2023).

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 (${\rm\Lambda 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
Ade  _et al._ (2016). 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 Hui and Others
(2017). Cold-DM-only simulation showed that the DM density has a cusp $\propto
1/r$ in the halo center, which is a known feature of the so-called Navarro-
Frenk-White (NFW) profile Navarro  _et al._ (1997). 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 Hui
and Others (2017). Observations of the galactic rotation curve have shown both
evidences for cusp and core-like distributions Moore  _et al._ (1999). 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 Heißel  _et al._ (2022); Lacroix (2018);
Shen  _et al._ (2024); Zakharov  _et al._ (2007). 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
Chakrabarti  _et al._ (2021). 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 Dubinski and Carlberg
(1991).

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 Damour and Deruelle
(1985). We invoke PN equation of motion for the two-body orbital dynamics
Barker and O’Connell (1975); Kidder (1995); Blanchet (2014),

$\ddot{\bm{r}}\equiv\frac{{\rm d}^{2}\bm{r}}{{\rm d}t^{2}}=\ddot{\bm{r}}_{\rm
N}+\ddot{\bm{r}}_{\rm 1PN}+\ddot{\bm{r}}_{\rm SO}+\ddot{\bm{r}}_{\rm
Q}+\ddot{\bm{r}}_{\rm 2PN}+\ddot{\bm{r}}_{\rm 2.5PN}+\cdots\,,$ (2)

where $\bm{r}\equiv\bm{r}_{\rm PSR}-\bm{r}_{\bullet}$ 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 $\ddot{\bm{r}}_{\rm N}=-GM\hat{\bm{n}}/r^{2}$, with $M\equiv
M_{\bullet}+m_{\rm PSR}$, $r\equiv\left|\bm{r}\right|$, and
$\hat{\bm{n}}\equiv\bm{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,

$\displaystyle\ddot{\bm{r}}_{\rm 1PN}$ $\displaystyle=$
$\displaystyle-\frac{GM}{c^{2}r^{2}}\left\\{\left[\left(-4-2\eta\right)\frac{GM}{r}+\left(1+3\eta\right)v^{2}-\frac{3}{2}\eta\dot{r}^{2}\right]\hat{\bm{n}}-\left(4-2\eta\right)\dot{r}\bm{v}\right\\}\,,$
(3a) $\displaystyle\ddot{\bm{r}}_{\rm SO}$ $\displaystyle=$
$\displaystyle\chi_{\bullet}\frac{G^{2}M^{2}}{4c^{3}r^{3}}\left(1+\sqrt{1-4\eta}\right)\left\\{\left[12\hat{\bm{s}}\cdot\left(\hat{\bm{n}}\times\bm{v}\right)\right]\hat{\bm{n}}+\left[\left(9+3\sqrt{1-4\eta}\right)\dot{r}\right]\left(\hat{\bm{n}}\times\hat{\bm{s}}\right)-\left(7+\sqrt{1-4\eta}\right)\left(\bm{v}\times\hat{\bm{s}}\right)\right\\}\,,$
(3b) $\displaystyle\ddot{\bm{r}}_{\rm Q}$ $\displaystyle=$ $\displaystyle-
q_{\bullet}\frac{3G^{3}M_{\bullet}^{2}M}{2c^{4}r^{4}}\left\\{\left[5\left(\hat{\bm{n}}\cdot\hat{\bm{s}}\right)^{2}-1\right]\hat{\bm{n}}-2\left(\hat{\bm{n}}\cdot\hat{\bm{s}}\right)\hat{\bm{s}}\right\\}\,,$
(3c)

where $\dot{r}\equiv{\rm d}r/{\rm d}t$, $\bm{v}\equiv{\rm d}\bm{r}/{\rm d}t$,
$v\equiv\left|\bm{v}\right|$, $\eta\equiv m_{\rm PSR}M_{\bullet}/M^{2}$, and
$\hat{\bm{s}}$ is the unit vector pointing along the BH spin. In particular,
in our simulation the acceleration at 2PN order, $\ddot{\bm{r}}_{\rm 2PN}$,
which is smaller by ${\cal O}\left(v^{2}/c^{2}\right)$ than
$\ddot{\bm{r}}_{\rm 1PN}$, and the acceleration at 2.5PN order, which leads to
the emission of gravitational waves, are not included Blanchet (2014). The
acceleration caused by the BH’s quadrupole is numerically at the 2PN order.
For consistency we should include at least the $\ddot{\bm{r}}_{\rm 2PN}$ 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_{\rm PSR}/M_{\bullet}<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 $\eta$ Barker and O’Connell (1975); 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,
$\dot{\bm{r}}_{\bullet}=0$ and $\bm{v}=\bm{v}_{\rm PSR}$. With $\eta\simeq 0$,
Eq. (3) simplifies to,

$\displaystyle\ddot{\bm{r}}_{\rm 1PN}$ $\displaystyle\simeq$
$\displaystyle-\frac{GM_{\bullet}}{c^{2}r^{2}}\left[\left(-\frac{4GM}{r}+v^{2}\right)\hat{\bm{n}}-4\dot{r}\bm{v}\right]\,,$
(4a) $\displaystyle\ddot{\bm{r}}_{\rm SO}$ $\displaystyle\simeq$
$\displaystyle\chi_{\bullet}\frac{6G^{2}M_{\bullet}^{2}}{c^{3}r^{3}}\left[\hat{\bm{s}}\cdot\left(\hat{\bm{n}}\times\bm{v}\right)\hat{\bm{n}}+\dot{r}\left(\hat{\bm{n}}\times\hat{\bm{s}}\right)-\frac{2}{3}\left(\bm{v}\times\hat{\bm{s}}\right)\right]\,,$
$\displaystyle\ddot{\bm{r}}_{\rm Q}$ $\displaystyle\simeq$ $\displaystyle-
q_{\bullet}\frac{3G^{3}M_{\bullet}^{3}}{2c^{4}r^{4}}\left\\{\left[5\left(\hat{\bm{n}}\cdot\hat{\bm{s}}\right)^{2}-1\right]\hat{\bm{n}}-2\left(\hat{\bm{n}}\cdot\hat{\bm{s}}\right)\hat{\bm{s}}\right\\}\,.$
(4c)

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.

Figure 1:  Coordinate system and notations that are used for a pulsar orbit
around Sgr A*. The frame
$\left(\hat{\bm{I}}_{0},\hat{\bm{J}}_{0},\hat{\bm{K}}_{0}\right)$ has
$\hat{\bm{K}}_{0}$ pointing from the Earth towards Sgr A*, and
$\left(\hat{\bm{I}}_{0},\hat{\bm{J}}_{0}\right)$ constitutes the sky plane
with $\hat{\bm{I}}_{0}$ pointing towards east and $\hat{\bm{J}}_{0}$ towards
north. The orientation of the pulsar orbit is determined by the inclination,
$i$, the longitude of periastron, $\omega$, and the longitude of ascending
node, $\Omega$. The Sgr A* sits at the origin, and its spin points towards
$\hat{\bm{s}}\equiv\left(\sin\lambda_{\bullet}\cos\eta_{\bullet},\,\sin\lambda_{\bullet}\sin\eta_{\bullet},\,\cos\lambda_{\bullet}\right)$
in the coordinate system; in the figure $\lambda_{\bullet}$ and
$\eta_{\bullet}$ are not shown for clarity. The position of the pulsar at the
reference time $t=0$ is $\bm{r}_{0}$, and it is described by the true anomaly
$\theta_{0}$.

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, $\bm{r}$, and instantaneous velocity, $\bm{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,

$\displaystyle M_{\bullet}$ $\displaystyle=$ $\displaystyle 4.3\times
10^{6}\,{\rm M}_{\odot}\,,$ (5a) $\displaystyle\chi_{\bullet}$
$\displaystyle=$ $\displaystyle
0.6\,,\quad\lambda_{\bullet}=\frac{1}{3}\pi\,,\quad\eta_{\bullet}=\frac{5}{9}\pi\,,\quad
q_{\bullet}=-0.36\,,$ (5b) $\displaystyle P_{b}$ $\displaystyle=$
$\displaystyle 0.5\,{\rm yr}\,,\quad e=0.8\,,$ (5c) $\displaystyle i$
$\displaystyle=$
$\displaystyle\frac{1}{5}\pi\,,\quad\omega=\frac{5}{7}\pi\,,\quad\theta_{0}=\frac{1}{3}\pi\,.$
(5d)

The angles are quasi-randomly chosen, to avoid special orientations that might
render the parameter-estimation problem degenerate.111The value of $\lambda$
was inspired by Psaltis  _et al._ (2015). We integrate the orbit for $t_{\rm
total}=5\,$yr, that contains $t_{\rm total}/P_{b}\simeq 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
$\Omega=0$. This choice is generic, because in this case the rotation of the
pulsar orbit around $\hat{\bm{K}}_{0}$ does not affect TOAs if we ignore the
proper motion of the Sgr A* Taylor (1994).

Figure 2: Illustration of the pulsar orbits around Sgr A* with $P_{b}=0.5$ yr,
$e=0.8$, and $\chi_{\bullet}=0.6$; see text and Eq. (5) for more details.

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 $\ddot{\bm{r}}_{\rm
1PN}$ and the frame-dragging effects from the BH spin Lense and Thirring
(1918); Schiff (1960).

##  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 Blandford and Teukolsky
(1976); Damour and Deruelle (1986); Taylor (1994); Shao (2023); Hu  _et al._
(2023). 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,

$N(T)=N_{0}+\nu T+\frac{1}{2}\dot{\nu}T^{2}\,,$ (6)

where $\nu\equiv 1/P$ is the pulsar’s spin frequency, $\dot{\nu}$ is the time
derivative of $\nu$, and $T$ is the proper time of the pulsar Damour and
Deruelle (1986). 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,

$\frac{{\rm d}T}{{\rm
d}t}=1-\frac{GM_{\bullet}}{c^{2}r}-\frac{v^{2}}{2c^{2}}+\cdots\,.$ (7)

Integrating the above equation gives an Einstein delay that accounts for the
gravitational redshift and special-relativistic time-dilation effects
Blandford and Teukolsky (1976); Damour and Deruelle (1986),

$\Delta_{\rm E}\equiv t-T\,.$ (8)

To drop a term that is linear in time, one can redefine the pulsar spin Damour
and Deruelle (1986) and Eq. (7) takes a form as

$\frac{{\rm d}T}{{\rm
d}t}=\frac{1-\frac{GM_{\bullet}}{c^{2}r}-\frac{v^{2}}{2c^{2}}}{\left\langle
1-\frac{GM_{\bullet}}{c^{2}r}-\frac{v^{2}}{2c^{2}}\right\rangle}\,.$ (9)

Roughly speaking, the $\left\langle\,\cdot\right\rangle$ term that appears in
the denominator means to average over the pulsar orbital period, so that
$\left\langle{\rm d}T/{\rm d}t\right\rangle=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, $\Delta_{\rm E}$ is found to be Blandford and Teukolsky
(1976); Damour and Deruelle (1986),

$\Delta_{\rm
E}=\frac{2e}{c^{2}}\left(\frac{G^{2}M_{\bullet}^{2}P_{b}}{2\pi}\right)^{1/3}\sin
u\,,$ (10)

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,

$\Delta_{\rm R}\equiv\frac{1}{c}\hat{\bm{K}}_{0}\cdot\bm{r}\,.$ (11)

We also include the lowest-order propagation time delay caused by the
curvature of the Sgr A*, called the 1PN Shapiro delay Shapiro (1964);
Blandford and Teukolsky (1976),

$\Delta_{\rm
S}=-\frac{2GM_{\bullet}}{c^{3}}\log\left(r-\bm{r}\cdot\hat{\bm{K}}_{0}\right)\,.$
(12)

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 Damour and Deruelle (1986),

$t^{\rm TOA}=T+\Delta_{\rm R}+\Delta_{\rm E}+\Delta_{\rm S}\,.$ (13)

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 (1997) and Wex and Kopeikin (1999) 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.

Figure 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
$\Delta_{\rm E}$ is removed; it can be absorbed into a redefinition of
pulsar’s spin.

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 Manchester  _et
al._ (2005); Liu  _et al._ (2012). Notice that for our fiducial orbit, the
inclination is only $\pi/5=36^{\circ}$. Even with such a small inclination,
the Shapiro delay is already numerous; this is also true even for face-on
orbits Liu  _et al._ (2012).

##  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 $\bm{\Theta}$ instead of the
pulsar timing model described in Sec. III that calculates $t^{\rm
TOA}(N;\bm{\Theta})$ Damour and Deruelle (1986). In principle, giving the
system’s parameters $\bm{\Theta}$ and the observed TOA, $t^{\rm TOA}$, we can
first integrate the pulsar’s orbital motion and get a series of time delays
that are related to the coordinate time: $\Delta_{\rm R}(t_{i})\,,\Delta_{\rm
S}(t_{i})\,,\Delta_{\rm E}(t_{i})$. Then we can do interpolation for these
series and solve the implicit equation, $t^{\rm TOA}=t+\Delta_{\rm
R}(t)+\Delta_{\rm S}(t)$, to get the coordinate pulse emission time, $t$.
Finally we use $T=t-\Delta_{\rm 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^{\rm TOA}$ that is directly related to the
observation.

From Eq. (8) and Eq. (13), we have

$t^{\rm TOA}=t+\Delta_{\rm R}(t)+\Delta_{\rm S}(t)\,.$ (14)

Taking a derivative on both sides and insert the explicit forms shown in Eq.
(11) and Eq. (12), we get

$\frac{{\rm d}t^{\rm TOA}}{{\rm
d}t}=1+\frac{1}{c}\hat{\bm{K}}_{0}\cdot\bm{v}-\frac{2GM_{\bullet}}{c^{3}}\frac{\hat{\bm{n}}\cdot\bm{v}-\bm{v}\cdot\hat{\bm{K}}_{0}}{r-\bm{r}\cdot\hat{\bm{K}}_{0}}\,.$
(15)

Combine this with Eq. (3) and Eq. (7), we can get the complete differential
equations used in the inverse timing formula, which are

$\displaystyle\frac{{\rm d}\bm{r}}{{\rm d}t^{\rm TOA}}$ $\displaystyle=$
$\displaystyle\bm{v}\frac{{\rm d}t}{{\rm d}t^{\rm TOA}}\,,$ (16)
$\displaystyle\frac{{\rm d}\bm{v}}{{\rm d}t^{\rm TOA}}$ $\displaystyle=$
$\displaystyle\ddot{\bm{r}}\frac{{\rm d}t}{{\rm d}t^{\rm TOA}}\,,$ (17)
$\displaystyle\frac{{\rm d}t}{{\rm d}t^{\rm TOA}}$ $\displaystyle=$
$\displaystyle\frac{{\rm d}t}{{\rm d}t^{\rm TOA}}\,,$ (18)
$\displaystyle\frac{{\rm d}\Delta_{\rm E}}{{\rm d}t^{\rm TOA}}$
$\displaystyle=$ $\displaystyle\left(1-\frac{{\rm d}T}{{\rm
d}t}\right)\frac{{\rm d}t}{{\rm d}t^{\rm TOA}}\,,$ (19)

where ${\rm d}t/{\rm d}t^{\rm 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 $\bm{\Theta}$ and $\Delta_{\rm E}=0$, we can
simply calculate $t^{\rm TOA}_{0}=\Delta_{\rm R}\big{|}_{t=0}+\Delta_{\rm
S}\big{|}_{t=0}$. The initial conditions now are at $t^{\rm TOA}=t^{\rm
TOA}_{0}$, and the system has parameters $\bm{\Theta}$ and $\Delta_{\rm E}=0$.

Integrating the above equations from $t^{\rm TOA}_{0}$ to $t^{\rm TOA}$, we
can get the related pulsar proper time $T$ by

$T=t^{\rm TOA}-\Delta_{\rm R}\big{(}t^{\rm TOA}\big{)}-\Delta_{\rm
E}\big{(}t^{\rm TOA}\big{)}-\Delta_{\rm S}\big{(}t^{\rm TOA}\big{)}\,,$ (20)

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 $\bm{\Theta}$, are given, one obtains
the pulsar rotation number $N$ as a function of $t^{\rm TOA}$ without any
ambigiuity. The parameters include,

$\bm{\Theta}=\bm{\Theta}_{\bullet}\cup\bm{\Theta}_{\rm
orbit}\cup\bm{\Theta}_{\rm PSR}\cup\bm{\Theta}_{\rm pert}\,,$ (21)

where $\bm{\Theta}_{\bullet}$ includes parameters of the BH in Eqs. (5a–5b),
$\bm{\Theta}_{\rm orbit}$ includes parameters of the orbit in Eqs. (5c–5d),
$\bm{\Theta}_{\rm PSR}$ includes parameters of the pulsar spin in Eq. (6), and
$\bm{\Theta}_{\rm pert}$ includes parameters from the perturbation. Therefore,
we denote the pulsar rotation number as $N\left(\bm{\Theta};t^{\rm
TOA}\right)$. $\bm{\Theta}_{\rm pert}$ is to be introduced in Sec. VI, and we
omit it in this section.

Here in $\bm{\Theta}_{\bullet}$, we treat $\chi_{\bullet}$ and $q_{\bullet}$
as independent variables. In GR, the no-hair theorem poses
$q_{\bullet}=-\chi_{\bullet}^{2}$ Israel (1967); Chrusciel  _et al._ (2012).
Our treatment can be viewed as an expansion of the multipoles of the spacetime
Thorne (1980) with the three lowest-order moments, namely the monopole
$M_{\bullet}$, the dipole $S=\chi_{\bullet}GM^{2}/c$, and the quadrupole
$Q=q_{\bullet}G^{2}M^{3}/c^{4}$. The post-Newtonian expansion in Eq. (3)
provides the possibility to treat $\chi_{\bullet}$ and $q_{\bullet}$
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 $\bm{\Theta}$ is,

$P\left(\bm{\Theta}\left|t^{\rm
TOA}\right.\right)\propto\exp\left(-\frac{P^{2}}{2}\sum_{i=1}^{N_{\rm
TOA}}\frac{\left[N^{(i)}\left(\bm{\Theta}\right)-N^{(i)}\left(\bar{\bm{\Theta}}\right)\right]^{2}}{\sigma_{\rm
TOA}^{2}}\right)\,,$ (22)

where $N^{(i)}\left(\bm{\Theta}\right)\equiv N\left(\bm{\Theta};t^{\rm
TOA}_{i}\right)$, 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. Shao (2016,
2017). Nevertheless, it is statistically equivalent to use the noiseless
templates, $N\left(\bar{\bm{\Theta}};t^{\rm TOA}\right)$ in Eq. (22). The
addition of random noises only shifts $\ln P\left(\bm{\Theta}\left|t^{\rm
TOA}\right.\right)$ 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 (2017) for simulations of
parameter estimation with gravitational-wave waveform templates.

As a standard procedure for parameter estimation Damour and Deruelle (1986);
Edwards  _et al._ (2006), we estimate the measurement uncertainties of these
parameters via the covariance matrix

$C_{\alpha\beta}=\left(\frac{\partial^{2}\mathcal{L}}{\partial\Theta^{\alpha}\partial\Theta^{\beta}}\right)^{-1}\,,$
(23)

where $\mathcal{L}=-\ln P\big{(}\Theta|t^{\rm TOA}\big{)}$ is the log-
likelihood function. We assume the timing precision $\sigma_{\rm TOA}=1\,{\rm
ms}$ in our simulation, which is a relatively conservative estimate for future
observations Liu  _et al._ (2012).

Figure 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 $\nu=2\,{\rm Hz}$ and the timing precision
$\sigma_{\rm TOA}=10^{-3}\,{\rm s}$. Simulations include weekly observations
over a 5-yr interval.

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._ (2012). Previous
studies suggest that $\sim 10^{3}$ pulsars can be expected to be orbiting
around Sgr A* with $P_{b}<100\,{\rm yr}$ Wharton  _et al._ (2012), while the
innermost some of them could be in orbits as tight as $\sim 100-500\,{\rm au}$
from Sgr A* Zhang  _et al._ (2014), which correspond to the orbital periods
considered here. From Fig. 4, we can see that for orbits with $P_{b}\lesssim
0.5\,{\rm yr}$, the measurement of the spin and quadrupole parameters can be
better than 1%, which is consistent with the results in Refs. Liu  _et al._
(2012); Bower  _et al._ (2018).

Figure 5: 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.

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
$\big{(}1-e^{2}\big{)}$ in the logarithmic scale. This can be explained
intuitively by that the secular effects of mass, spin, and quadrupole are
proportional to $\big{(}1-e^{2}\big{)}^{-1}$, $\big{(}1-e^{2}\big{)}^{-3/2}$,
and $\big{(}1-e^{2}\big{)}^{-2}$ respectively Wex and Kopeikin (1999); Barker
and O’Connell (1975). From this figure we can also conclude that for obits
with very high eccentricities, say $e\gtrsim 0.8$, the precision of the
quadrupole determination can reach the precision of the spin determination or
even better, which comes from the $\propto r^{-4}$ nature of the quadrupole
interaction.

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

As shown in Refs. Liu  _et al._ (2012); Zhang and Saha (2017), there are
leading-order degeneracies among spin parameters. The observable secular
effects caused by the BH spin including the advance of the periastron $\omega$
and the change of the inclination angle $i$ Wex and Kopeikin (1999). 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 $\ddot{\omega}$, or periodic effects are needed when
measuring the BH spin and its orientation Liu  _et al._ (2012). In Fig. 6 we
plot the correlation between the spin parameters, where blue lines give the
predicted leading order degeneracies determined by Zhang and Saha (2017)

$\displaystyle\mathbb{S}\cdot\big{(}\hat{\bm{L}}\times\hat{\bm{K}}_{0}\big{)}$
$\displaystyle=$ $\displaystyle{\rm const.}\,,$ (24) $\displaystyle
2\big{(}\hat{\bm{L}}\cdot\mathbb{S}\big{)}+\big{(}\hat{\bm{A}}\cdot\mathbb{S}\big{)}\frac{\cos
i}{\sin\omega\sin i}$ $\displaystyle=$ $\displaystyle{\rm const.}\,,$ (25)

where $\mathbb{S}=S\hat{\bm{s}}$ is the spin vector of the BH and
$\hat{\bm{A}}$ 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 Merritt  _et al._ (2010). 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 Psaltis  _et al._ (2016).
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.1 and Sec. VI.2. In Sec. VI.3 we present the extended timing model and the
parameter estimation results.

###  VI.1 The Generalized NFW Profile

The standard $\rm\Lambda 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 Navarro
_et al._ (1997)

$\rho_{\rm NFW}(r)=\frac{\rho_{0}}{(r/R_{s})(1+r/R_{s})^{2}}\,,$ (26)

where $R_{s}$ is the scale radius and $\rho_{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}\simeq 20\,{\rm kpc}$, and $\rho_{0}\big{|}_{r\approx 8\,{\rm
kpc}}=0.4\,{\rm GeV\,cm^{-3}}$ at the location of our Solar System McMillan
(2017).

A feature of the NFW profile is the singular density profile at the center
where $\rho(r)\propto r^{-1}$, which is called a density cusp. However,
observations of the rotation curve of low-surface-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 Moore  _et al._ (1999);
Vargya  _et al._ (2022). 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 Hui and Others (2017). 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 Pontzen and Governato (2012); Cole  _et al._ (2011), 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 Vargya  _et al._ (2022). Nevertheless, to account for
different model predictions, here we consider a generalized model, the gNFW
model Diemand  _et al._ (2008)

$\rho_{\rm
gNFW}(r)=\frac{\rho_{0}}{(r/R_{s})^{\gamma}(1+r/R_{s})^{3-\gamma}}\,.$ (27)

The central behaviors of different models are characterized by the power-law
index $\gamma$.

Gondolo and Silk (1999) 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}\equiv GM_{\bullet}/v_{0}^{2}$, which is
about $1.7\,{\rm pc}$ for Sgr A* Gultekin  _et al._ (2009). For collisionless
DM with a polytrope phase-space distribution, a DM spike with $\rho_{\rm
sp}(r)\propto r^{-\gamma_{\rm sp}}$ will form, with $\gamma_{\rm
sp}=(9-2\gamma)/(4-\gamma)$. For a pulsar with orbital period $P_{b}\sim
1\,{\rm yr}$ has a semi-major axis $\sim 1\,{\rm 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 Gondolo
and Silk (1999). If so, due to the very high DM densities in the halo center,
it may act as a strong gamma-ray source Gondolo and Silk (1999), 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 $\rho_{\rm in}(r)\propto
r^{-\gamma_{\rm in}}$, where $\gamma_{\rm in}\simeq 0.5$ Vasiliev (2007). Such
a weak annihilation cusp will form around a radius $R_{\rm in}$ where the DM
density reaches a critical density $\rho_{\rm ann}$.

Combining the above discussion, we consider the DM density profile as below
Shao  _et al._ (2018),

$\rho_{\rm DM}(r)=\left\\{\begin{aligned} &\frac{\rho_{\rm sp}(r)\rho_{\rm
in}(r)}{\rho_{\rm sp}(r)+\rho_{\rm in}(r)}\,,&4GM_{\bullet}/c^{2}\leq r<R_{\rm
sp}\\\ &\rho_{\rm gNFW}(r)\,,&r\geq R_{\rm sp}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \end{aligned}\right.\,.$ (28)

If there is no DM annihilation, the profile is

$\rho_{\rm DM}(r)=\left\\{\begin{aligned} &\rho_{\rm
sp}(r)&4GM_{\bullet}/c^{2}\leq r<R_{\rm sp}\\\ &\rho_{\rm gNFW}(r)\,,&r\geq
R_{\rm sp}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \end{aligned}\right.\,,$ (29)

where $R_{\rm sp}$ is chosen to equal to the gravitational radius $R_{h}$.
Densities $\rho_{\rm sp}(r)$ and $\rho_{\rm in}(r)$ are determined by the
continuous condition at $R_{\rm sp}$. For the model with annihilation, we
assume $\rho_{\rm sp}(R_{\rm in})=\rho_{\rm in}(R_{\rm in})=2\rho_{\rm ann}$
to interpolate these two profiles.

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

###  VI.2 The Einasto Profile

The Einasto profile is also a commonly used DM distribution in which the
logarithmic density slope shows a power-law behavior Einasto (1965); Navarro
_et al._ (2004); Merritt  _et al._ (2005), 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 Navarro  _et al._ (2004); Wang  _et al._ (2020). The DM
density profile of this model can be written as Shen  _et al._ (2024)

$\rho_{\rm
E}(r)=\rho_{0}\exp\left\\{-\frac{2}{\alpha}\left[\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right]\right\\}\,,$
(30)

where $\rho_{0}$ is the DM density at the scale radius $r_{s}$ and $\alpha$ is
the inverse of the Einasto index. Fitting the rotation curve and globular
cluster kinematics from Gaia data gives Wang  _et al._ (2022)

$\displaystyle\rho_{\rm E}(8.2\,{\rm kpc})$
$\displaystyle=0.008\,M_{\odot}\,{\rm pc^{-3}}\,,$ (31) $\displaystyle r_{s}$
$\displaystyle=12\,{\rm kpc}\,,$ (32) $\displaystyle\alpha$
$\displaystyle=0.32\,.$ (33)

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 circular-orbit approximation Shen  _et al._ (2024). Following
Shen  _et al._ (2024), we estimate the DM spike via

$\displaystyle r_{i}M_{{\rm tot},\,i}(r_{i})$ $\displaystyle=r_{f}M_{{\rm
tot},\,f}(r_{f})\,,$ (34) $\displaystyle M_{{\rm DM},\,i}(r_{i})$
$\displaystyle=M_{{\rm DM},\,f}(r_{f})\,,$ (35)

where $M_{\rm DM}(r)$ is the mass of DM enclosed in the radius $r$

$M_{\rm DM}(r)=\int_{0}^{r}4\pi r^{2}\rho_{\rm DM}(r){\rm d}r\,,$ (36)

and $M_{\rm tot}(r)=M_{\rm SMBH}+M_{\rm DM}(r)$. Exact distribution of the DM
spike in the Einasto model can be numerically calculated as shown in Ref. Shen
_et al._ (2024). 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
$\rho_{E}(r)=\rho_{0}e^{2/\alpha}$. Combining with above equations, one
derives that the DM spike density in the Einasto model is $\rho_{\rm
E,\,sp}(r)\propto r^{9/4}$, which is effectively to have $\gamma=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
Shen  _et al._ (2024). 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
Poisson and Will (2014), 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.

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

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_{\odot}$, while for the DM model
with annihilation, timing a pulsar with an orbital period $P_{b}\lesssim
10\,{\rm yr}$ can only give the information of the weak cusp caused by the DM
annihilation, described by $\gamma_{\rm 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 $\gamma_{\rm 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 $\gamma_{\rm sp}$, which is related to the
underlying DM profile via $\gamma_{\rm sp}=(9-2\gamma)/(4-\gamma)$.

###  VI.3 Prameter Estimation

Compared to the mass of the central BH, the DM mass is still very small (about
$10^{-5}\,M_{\bullet}$ 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.

Figure 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). Figure 10: Timing residuals caused by the DM perturbation before
fitting. We use the DM model with spike but no annihilation, and we set
$R_{\rm sp}=R_{h}$. The system parameters are given in Eqs. (5a–5d). Figure
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.

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
$\bar{\bm{\Theta}}$ 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_{\odot}$ Abuter  _et al._ (2022); Heißel  _et al._ (2022), 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 $\lesssim 1000\,M_{\odot}$, although the star clusters will
also contribute to the extended mass in this scale Heißel  _et al._ (2022).

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

$\rho_{\rm
DM}(r)=\rho_{0}\left(\frac{r}{4GM_{\bullet}/c^{2}}\right)^{-\gamma_{\rm
sp}}\,,\ \ r\geq 4GM_{\bullet}/c^{2}\,,$ (37)

where $\rho_{0}$ now represents the central density of the DM spike and
$\gamma_{\rm sp}$ is related to the original DM profile, where we use the NFW
profile with $\gamma=1$ in our simulation.

Figure 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 $\rho_{0}$ and $\gamma_{\rm
sp}$ are the parameters for the NFW model with spike.

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}\gtrsim 0.5\,{\rm
yr}$ and an orbital eccentricity $e\sim 0.8$ can give a fractional measurement
uncertainty of $\sim 1\%$ in $\gamma_{\rm sp}$ for the spike model, which is
related to a $\sim 20\%$ fractional precision in $\gamma$. This result is
comparable to the constraints from fitting kinematic data of maser and other
observations, which give $\gamma=0.79\pm 0.32$ for our Galaxy McMillan (2017).
However, by timing a pulsar around Sgr A*, we can constrain the DM structure
in the length scale of $\sim{\rm 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 $\sim{\rm kpc}$ scales.

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

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\to 0$. Only when the pulsar moves at different
radii it can sense the DM’s radial distribution.

In the above, we have used unrelated $\chi_{\bullet}$ and $q_{\bullet}$,
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_{\bullet}=-\chi_{\bullet}^{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._
(2022) that, the orbital features caused by the extended mass distribution are
significant in the orbital section $\theta\in(5\pi/6,\,7\pi/6)$, which is
different from the spin or quadrupole effects that are most significant near
the periastron. So setting $q_{\bullet}=-\chi_{\bullet}^{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 leading-
order effects of various time delays being taken into account. Our simulations
show that for a pulsar with an orbital period $P_{b}\lesssim 0.5\,{\rm yr}$
and an orbital eccentricity $e\sim 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 Liu  _et
al._ (2012); Psaltis  _et al._ (2016) or a fully general-relativistic
treatment Zhang and Saha (2017). 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}\gtrsim 0.5\,{\rm yr}$, and an orbital
eccentricity $e\sim 0.8$ can provide a 1% measurement of the power-law index
$\gamma_{\rm sp}$ of the DM spike, which relates to a measurement of the
underlying DM model with a $~{}20\%$ precision in the index $\gamma$. Such a
precision is comparable with, but complementary to, the observations at the
Galactic scale, usually done at $\sim$ kpc scales, which can be an equilibrium
based kinematic analysis McMillan (2017), or some extreme-precision time-
series measurements of Galactic accelerations Chakrabarti  _et al._ (2020,
2022). 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.

###### Acknowledgements.

We thank Norbert Wex for helpful discussions, and Kuo Liu for carefully
reading the manuscript. This work was supported by the National SKA Program of
China (2020SKA0120300), the National Natural Science Foundation of China
(11991053, 11975027, 12273006), the Max Planck Partner Group Program funded by
the Max Planck Society, and the High-Performance Computing Platform of Peking
University.

## References

  * Hawking and Ellis (2023) S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2023).
  * Chandrasekhar (1983) S. Chandrasekhar,  _The Mathematical Theory of Black Holes_ (Oxford University Press, 1983).
  * Chrusciel  _et al._ (2012) P. T. Chrusciel, J. Lopes Costa, and M. Heusler, Living Rev. Rel. 15, 7 (2012).
  * Unruh and Wald (2017) W. G. Unruh and R. M. Wald, Rept. Prog. Phys. 80, 092002 (2017).
  * McConnell and Ma (2013) N. J. McConnell and C.-P. Ma, Astrophys. J. 764, 184 (2013).
  * Wex and Kopeikin (1999) N. Wex and S. Kopeikin, Astrophys. J. 514, 388 (1999).
  * Will (2008) C. M. Will, Astrophys. J. Lett. 674, L25 (2008).
  * Ghez  _et al._ (2008) A. M. Ghez _et al._ , Astrophys. J. 689, 1044 (2008).
  * Genzel  _et al._ (2010) R. Genzel, F. Eisenhauer, and S. Gillessen, Rev. Mod. Phys. 82, 3121 (2010).
  * Akiyama  _et al._ (2022a) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 930, L12 (2022a).
  * Akiyama  _et al._ (2022b) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 930, L17 (2022b).
  * Israel (1967) W. Israel, Phys. Rev. 164, 1776 (1967).
  * Carter (1971) B. Carter, Phys. Rev. Lett. 26, 331 (1971).
  * Robinson (1975) D. C. Robinson, Phys. Rev. Lett. 34, 905 (1975).
  * Hansen (1974) R. O. Hansen, J. Math. Phys. 15, 46 (1974).
  * Thorne (1980) K. S. Thorne, Rev. Mod. Phys. 52, 299 (1980).
  * Kramer  _et al._ (2004) M. Kramer, D. C. Backer, J. M. Cordes, T. J. W. Lazio, B. W. Stappers, and S. Johnston, New Astron. Rev. 48, 993 (2004), arXiv:astro-ph/0409379 .
  * Liu  _et al._ (2012) K. Liu _et al._ , Astrophys. J. 747, 1 (2012).
  * Shao  _et al._ (2015) L. Shao _et al._ , PoS AASKA14, 042 (2015), arXiv:1501.00058 [astro-ph.HE] .
  * Zhang and Saha (2017) F. Zhang and P. Saha, Astrophys. J. 849, 33 (2017).
  * Bower  _et al._ (2018) G. C. Bower _et al._ , ASP Conf. Ser. 517, 793 (2018), arXiv:1810.06623 [astro-ph.HE] .
  * Weltman  _et al._ (2020) A. Weltman _et al._ , Publ. Astron. Soc. Austral. 37, e002 (2020), arXiv:1810.02680 [astro-ph.CO] .
  * Eatough  _et al._ (2023) R. P. Eatough _et al._ , in  _15th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Astrophysics, and Relativistic Field Theories_ (2023) arXiv:2306.01496 [astro-ph.HE] .
  * Abuter  _et al._ (2020) R. Abuter _et al._ (GRAVITY), Astron. Astrophys. 636, L5 (2020).
  * Abuter  _et al._ (2022) R. Abuter _et al._ (GRAVITY), Astron. Astrophys. 657, L12 (2022).
  * Heißel  _et al._ (2022) G. Heißel _et al._ , Astron. Astrophys. 660, A13 (2022).
  * Merritt  _et al._ (2010) D. Merritt _et al._ , Phys. Rev. D 81, 062002 (2010).
  * Akiyama  _et al._ (2022c) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 930, L15 (2022c).
  * Xu  _et al._ (2023a) R. Xu, D. Liang, and L. Shao, Phys. Rev. D 107, 024011 (2023a), arXiv:2209.02209 [gr-qc] .
  * Xu  _et al._ (2023b) R. Xu, D. Liang, and L. Shao, Astrophys. J. 945, 148 (2023b), arXiv:2302.05671 [gr-qc] .
  * Psaltis  _et al._ (2016) D. Psaltis _et al._ , Astrophys. J. 818, 121 (2016).
  * Dong  _et al._ (2022) Y. Dong, L. Shao, Z. Hu, X. Miao, and Z. Wang, JCAP 11, 051 (2022), arXiv:2210.16130 [astro-ph.HE] .
  * Kramer  _et al._ (2000) M. Kramer, B. Klein, D. R. Lorimer, P. Mueller, A. Jessner, and R. Wielebinski, ASP Conf. Ser. 202, 37 (2000).
  * Eatough  _et al._ (2013) R. P. Eatough _et al._ , Nature 501, 391 (2013).
  * Wharton  _et al._ (2012) R. S. Wharton, S. Chatterjee, J. M. Cordes, J. S. Deneva, and T. J. W. Lazio, Astrophys. J. 753, 108 (2012).
  * Siemion  _et al._ (2013) A. Siemion _et al._ , in  _Neutron Stars and Pulsars: Challenges and Opportunities after 80 years_, Vol. 291, edited by J. van Leeuwen (2013) pp. 57–57.
  * Liu  _et al._ (2021) K. Liu _et al._ , Astrophys. J. 914, 30 (2021), arXiv:2104.08986 [astro-ph.HE] .
  * Torne  _et al._ (2023) P. Torne _et al._ (EHT), Astrophys. J. 959, 14 (2023), arXiv:2308.15381 [astro-ph.HE] .
  * Cordes and Lazio (2002) J. M. Cordes and T. J. W. Lazio, (2002), arXiv:astro-ph/0207156 .
  * Wex (1995) N. Wex, Class. Quant. Grav. 12, 983 (1995).
  * Damour and Deruelle (1986) T. Damour and N. Deruelle, Ann. Inst. Henri Poincaré Phys. Théor. 44, 263 (1986).
  * Straub  _et al._ (2023) O. Straub _et al._ (GRAVITY), Astron. Astrophys. 672, A63 (2023).
  * Ade  _et al._ (2016) P. A. R. Ade _et al._ (Planck), Astron. Astrophys. 594, A13 (2016).
  * Hui and Others (2017) L. Hui and Others, Phys. Rev. D 95, 043541 (2017).
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997).
  * Moore  _et al._ (1999) B. Moore _et al._ , Mon. Not. Roy. Astron. Soc. 310, 1147 (1999).
  * Lacroix (2018) T. Lacroix, Astron. Astrophys. 619, A46 (2018), arXiv:1801.01308 [astro-ph.GA] .
  * Shen  _et al._ (2024) Z.-Q. Shen, G.-W. Yuan, C.-Z. Jiang, Y.-L. S. Tsai, Q. Yuan, and Y.-Z. Fan, Mon. Not. Roy. Astron. Soc. 527, 3196 (2024), arXiv:2303.09284 [astro-ph.GA] .
  * Zakharov  _et al._ (2007) A. F. Zakharov, A. A. Nucita, F. De Paolis, and G. Ingrosso, Phys. Rev. D 76, 062001 (2007).
  * Chakrabarti  _et al._ (2021) S. Chakrabarti _et al._ , Astrophys. J. Lett. 907, L26 (2021).
  * Dubinski and Carlberg (1991) J. Dubinski and R. G. Carlberg, Astrophys. J. 378, 496 (1991).
  * Damour and Deruelle (1985) T. Damour and N. Deruelle, Ann. Inst. Henri Poincaré Phys. Théor. 43, 107 (1985).
  * Barker and O’Connell (1975) B. M. Barker and R. F. O’Connell, Phys. Rev. D12, 329 (1975).
  * Kidder (1995) L. E. Kidder, Phys. Rev. D52, 821 (1995).
  * Blanchet (2014) L. Blanchet, Living Rev. Rel. 17, 2 (2014).
  * Psaltis  _et al._ (2015) D. Psaltis _et al._ , Astrophys. J. 798, 15 (2015).
  * Taylor (1994) J. H. Taylor, Rev. Mod. Phys. 66, 711 (1994).
  * Lense and Thirring (1918) J. Lense and H. Thirring, Phys. Z. 19, 156 (1918).
  * Schiff (1960) L. I. Schiff, Phys. Rev. Lett. 4, 215 (1960).
  * Blandford and Teukolsky (1976) R. Blandford and S. A. Teukolsky, Astrophys. J. 205, 580 (1976).
  * Shao (2023) L. Shao, Lect. Notes Phys. 1017, 385 (2023), arXiv:2206.15187 [gr-qc] .
  * Hu  _et al._ (2023) Z. Hu, X. Miao, and L. Shao, (2023), arXiv:2303.17185 [astro-ph.HE] .
  * Shapiro (1964) I. I. Shapiro, Phys. Rev. Lett. 13, 789 (1964).
  * Kopeikin (1997) S. M. Kopeikin, Journal of Mathematical Physics 38, 2587 (1997).
  * Manchester  _et al._ (2005) R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron. J. 129, 1993 (2005).
  * Shao (2016) L. Shao, Phys. Rev. D93, 084023 (2016).
  * Shao (2017) L. Shao, Class. Quant. Grav. 34, 175011 (2017).
  * Gaebel and Veitch (2017) S. M. Gaebel and J. Veitch, Class. Quant. Grav. 34, 174003 (2017).
  * Edwards  _et al._ (2006) R. T. Edwards, G. B. Hobbs, and R. N. Manchester, Mon. Not. Roy. Astron. Soc. 372, 1549 (2006).
  * Zhang  _et al._ (2014) F. Zhang, Y. Lu, and Q. Yu, Astrophys. J. 784, 106 (2014).
  * McMillan (2017) P. J. McMillan, Mon. Not. Roy. Astron. Soc. 465, 76 (2017).
  * Vargya  _et al._ (2022) D. Vargya _et al._ , Mon. Not. Roy. Astron. Soc. 516, 2389 (2022).
  * Pontzen and Governato (2012) A. Pontzen and F. Governato, Mon. Not. Roy. Astron. Soc. 421, 3464 (2012).
  * Cole  _et al._ (2011) D. Cole, W. Dehnen, and M. Wilkinson, Mon. Not. Roy. Astron. Soc. 416, 1118 (2011).
  * Diemand  _et al._ (2008) J. Diemand _et al._ , Nature 454, 735 (2008).
  * Gondolo and Silk (1999) P. Gondolo and J. Silk, Phys. Rev. Lett. 83, 1719 (1999).
  * Gultekin  _et al._ (2009) K. Gultekin _et al._ , Astrophys. J. 698, 198 (2009).
  * Vasiliev (2007) E. Vasiliev, Phys. Rev. D 76, 103532 (2007).
  * Shao  _et al._ (2018) L. Shao, N. Wex, and M. Kramer, Phys. Rev. Lett. 120, 241104 (2018).
  * Einasto (1965) J. Einasto, Trudy Astrofizicheskogo Instituta Alma-Ata 5, 87 (1965).
  * Navarro  _et al._ (2004) J. F. Navarro _et al._ , Mon. Not. Roy. Astron. Soc. 349, 1039 (2004).
  * Merritt  _et al._ (2005) D. Merritt _et al._ , Astrophys. J. Lett. 624, L85 (2005).
  * Wang  _et al._ (2020) J. Wang _et al._ , Nature 585, 39 (2020).
  * Wang  _et al._ (2022) J. Wang, F. Hammer, and Y. Yang, Mon. Not. Roy. Astron. Soc. 510, 2242 (2022).
  * Poisson and Will (2014) E. Poisson and C. M. Will,  _Gravity: Newtonian, Post-Newtonian, Relativistic_ (Cambridge University Press, 2014).
  * Chakrabarti  _et al._ (2020) S. Chakrabarti _et al._ , Astrophys. J. Lett. 902, L28 (2020).
  * Chakrabarti  _et al._ (2022) S. Chakrabarti _et al._ , Astrophys. J. Lett. 928, L17 (2022).