Exploring dark matter spike distribution around the Galactic centre with stellar orbits

Precise measurements of the stellar orbits around Sagittarius A* have established the existence of a supermassive black hole (SMBH) at the Galactic centre (GC). Due to the interplay between the SMBH and dark matter (DM), the DM density profile in the innermost region of the Galaxy, which is crucial for the DM indirect detection, is still an open question. Among the most popular models in the literature, the theoretical spike profile proposed by Gondolo and Silk (1999; GS hereafter) is well adopted. In this work, we investigate the DM spike profile using updated data from the Keck and VLT telescopes considering that the presence of such an extended mass component may affect the orbits of the S-stars in the Galactic center. We examine the radius and slope of the generalized NFW spike profile, analyze the Einasto spike, and discuss the influence of DM annihilation on the results. Our findings indicate that an initial slope of $\gamma \gtrsim 0.92$ for the generalized NFW spike profile is ruled out at a 95% confidence level. Additionally, the spike radius $R_{\rm sp}$ larger than 21.5 pc is rejected at 95% probability for the Einasto spike with $\alpha=0.17$, which also contradicts the GS spike model. The constraints with the VLT/GRAVITY upper limits are also projected. Although the GS NFW spike is well constrained by the Keck and VLT observation of S2, an NFW spike with a weak annihilation cusp may still be viable, as long as the DM annihilation cross section satisfies $\left<\sigma v \right>\gtrsim 7.7\times 10^{-27}~{\rm cm^3\,s^{-1}} (m_{\rm DM}/100~{\rm GeV})$ at 95% level.


INTRODUCTION
The nature of dark matter (DM), an invisible component that provides the additional gravitational force necessary to explain a range of phenomena across different scales, remains one of the most significant enigmas in the universe (Bergström 2000;Bertone & Hooper 2018).To date, the particle nature of DM remains largely unknown, despite numerous proposals for DM candidates in the scientific literature (Bertone et al. 2005;Feng 2010;Hu et al. 2000;de Laurentis et al. 2022).In addition, some systematic searches have been conducted to explore these candidates (Porter et al. 2011;Charles et al. 2016;Roszkowski et al. 2018;Liu et al. 2017;Zhu et al. 2022).Of particular interest in the indirect detection of DM is the Galactic centre (GC), where the density of DM peaks and certain excesses have been reported and investigated (e.g.Hooper & Goodenough 2011;Di Mauro 2021;Cholis et al. 2022;Bringmann et al. 2012;Zhou et al. 2015;Ackermann et al. 2015;Alemanno et al. 2022;Abe et al. 2023).However, it is challenging to determine the DM density in the inner Galaxy based on the rotation curve of interstellar gas, as outlined in previous studies (Sofue 2013).
The actual DM density profile may theoretically differ from the ★ E-mail: smingtsai@pmo.ac.cn (YLST), yzfan@pmo.ac.cn (YZF) DM-only halo as a result of the interplay between DM and the supermassive black hole (SMBH) (e.g.Detweiler 1980;Gondolo & Silk 1999;Yuan et al. 2022b;Cai et al. 2023).The DM spike, a DM structure even steeper than the cusp, can be formed along with the growth of the SMBH.The accumulation of DM particles in the GC can be comprehensively explained as a result of the inward flow of the ordinary interstellar medium, which undergoes accretion by the black hole (BH) due to frictional forces.This process leads to an increase in the mass and gravitational potential of the BH, ultimately leading to the accumulation of DM particles in the GC.According to the proposal by Gondolo & Silk (1999), if the SMBH in the GC grows adiabatically, the DM density could be enhanced by up to ten orders of magnitude.The profile (GS spike hereafter) is proportional to  − sp , with the spike slope range 2.25 <  sp < 2.5 (Gondolo & Silk 1999;Merritt et al. 2002;Sadeghian et al. 2013;Ferrer et al. 2017).This spike profile has been widely adopted in previous studies, leading to strong constraints on the DM annihilation rate (e.g.Gondolo 2000;Bertone et al. 2002;Fields et al. 2014;Lacroix et al. 2017;Xia et al. 2021;Liu et al. 2022;Balaji et al. 2023).
The formation of SMBHs is a complex process, and it may not always be adiabatic or initially located at the centre of the galaxy, leading to a flattening of the DM spike.Several studies have demonstrated that, by relaxing the ideal assumptions, the slope of the spike may differ.The slope of the spike can be weakened to 3/2 due to sufficient scatterings of DM particles with dense stellar populations (Gnedin & Primack 2004;Merritt 2004;Shapiro & Heggie 2022).The formation of SMBHs through an instantaneous gas collapse can result in a spike slope of 4/3 (Ullio et al. 2001).Major mergers of DM halos containing SMBHs can heat DM particles, producing a cusp with a slope of 0.5 within about 10 pc (Merritt et al. 2002).Additionally, the spike density can be weakened if the seed BH is massive enough and located off-centre (Ullio et al. 2001).However, other processes such as chaotic orbits in triaxial halos (Merritt & Poon 2004) or gravo-thermal collapse for self-interacting DM (Ostriker 2000) may enhance the spike.Still, the innermost history of the Galaxy is uncertain (Chen et al. 2023), making it difficult to determine the extent of influence on the DM spike.Therefore, in order to better constrain the DM parameters, it is necessary to probe the spike in observations.There are two options: the stellar orbits of S-stars (Weinberg et al. 2005;Lacroix 2018) and the gravitational wave (Eda et al. 2013;Li et al. 2022;Shen et al. 2023;Ghoshal & Strumia 2023), while our work only focus on the former one.
Over the past three decades, considerable efforts have been made to accurately measure the stellar kinematics in the innermost region of the Galaxy (Eckart & Genzel 1996;Ghez et al. 1998;Schödel et al. 2002;Ghez et al. 2003;Abuter et al. 2018a;Do et al. 2019).Thanks to the high resolution of the Keck observatory and the Very Large Telescope (VLT), the orbits of more than 40 S-stars are currently available (Gillessen et al. 2017;Peißker et al. 2020aPeißker et al. ,b, 2022;;Abuter et al. 2021a).These data have significantly improved our understanding of the SMBH properties (Genzel et al. 2010;Abuter et al. 2022;Zhang et al. 2015), the environment around the black hole (Amorim et al. 2019;Bar et al. 2019;Becerra-Vergara et al. 2021;Benisty & Davis 2022;Yuan et al. 2022a;Straub et al. 2023), and even the gravity theory (Hees et al. 2017;Abuter et al. 2020;Yan et al. 2022;Della Monica et al. 2023).Recently, due to the updated orbit measurements of S-stars (Gillessen et al. 2017;Do et al. 2019;Abuter et al. 2022) and the first image of Sgr A* presented by the Event Horizon Telescope (EHT) (Akiyama et al. 2022), interests in the DM spike are aroused.In a study by Lacroix (2018), the VLT and Keck observations of S2 based on the period before 2016 were analyzed, and the size of the spike was constrained.Their results only exclude the GS spike model for  ≳ 1.4 with 95% probability because the star had not yet reached its pericentre at that time.In a subsequent study by Abuter et al. (2020), the radius of the spike for  = 1 was further constrained by using the GRAVITY/VLT data.In another recent study by Nampalliwar et al. (2021), the Keck measurement of S2 was used to determine the inner radius and density of the spike, with a particular focus on the impact of the spike on the EHT shadow image.
In this work, we revisit and extend the analyses of DM spike with the public data of the S-stars from Keck and VLT.Our work not only updates the constraints on the radial extension parameter  sp for different values of , but also establishes, for the first time, limits on the spike slope  sp using the S2 orbit data.We further investigate the impact of combining data from multiple S-stars and calculate the Einasto profile as the initial density distribution. 1We also discuss the requirement of the annihilation cross section for the NFW spike with a weak cusp to escape the constraints of the S2 orbit.
This paper is organized as follows.In Sec. 2, we introduce the density profiles of DM spike.In Sec. 3, we present the orbit data of the S-stars, the post-Newtonian dynamical model, and the statistical method employed in our analysis.In Sec. 4, we show our constraints of the NFW spike and Einasto spike using the S2 orbital data from the Keck and VLT.In Sec. 5, we further discuss the constraints when more S-stars are combined, when a full-orbit GRAVITY measurements are adopted, and when DM particles can annihilate.Finally, we summarize our work in Sec. 6.

DARK MATTER PROFILE
It is a popular assumption that a spike profile can be formed after the adiabatic growth of SMBH (Gondolo & Silk 1999).Under such an adiabatic condition, one can obtain the analytical form of the NFW spike in the GS model assuming the conservation of angular momentum and radial action (Gondolo & Silk 1999;Young 1980).On the other hand, the analytical form of the Einasto spike in the GS model is hard to find, but one can obtain its circular-orbit approximation (Ullio et al. 2001;Blumenthal et al. 1986).In this section, we introduce the spike profiles: the NFW spike presented in Sec.2.1 and the Einasto spike presented in Sec.2.2.
Note that we are also phenomenologically interesting the general case of removing the assumptions in the GS model.Except presenting the GS spike profiles for both halo models, we will also compare generic spike with the GS spike by releasing the conditions and hence we can have one free parameter in the generic spike scenario.

The NFW spike
The most popular DM density model is the generalized Navarro-Frenk-White (gNFW) profile.The density at the Galactocentric radius  is (1) Here,  s is the scale radius with the condition d ln /d ln  = −2.The steepness of the density profile within  s is defined by , which equals 1 for the original NFW profile (Navarro et al. 1997).The density normalization  0 can be determined from the observations (McMillan 2017;Cautun et al. 2020;Benito et al. 2021;Wang et al. 2022).In our work, we adopt the gNFW profiles for different  from McMillan (2017) and mainly focus on the case of 0.5 ≤  ≤ 1.5.The upper value of  covers the results given in McMillan (2017) and Wang et al. (2022), and the lower value is adopted simply due to that the GS model with  < 0.5 is hardly constrained by the current data of the S-stars and is beyond our interests.
For the general gNFW spike profile (no adiabatic assumption), we adopt the piece-wise function (Lacroix 2018) where  sp0 =  gnfw,halo ( sp ) and  sch is the Schwarzschild radius of the GC SMBH which is 2  BH / 2 with  BH ≈ 4 × 10 6  ⊙ . sp and  sp are the spike radius and slope, respectively, which are two independent parameters to determine the spike density.
If considering the spike after the adiabatic growth of SMBH (GS spike), its spike slope is and the spike radius is where   is the scale factor interpolated from the values in Gondolo & Silk (1999).We can see that  GS sp and  GS sp both are the function of .As a reference, the predicted GS spike radii and slopes are given in Tab.B1.Upon comparison of the NFW halo profile with and without the spike profile, depicted by the red solid and dashed lines correspondingly in Fig. 1, a significant increase in the DM density of approximately six orders of magnitude can be observed at the pericentre of S2.

The Einasto spike
The Einasto profile is also a commonly used DM distribution model.By accounting for the power-law evolution of the logarithmic density slope with respect to the radius, this model can better fit the numeric simulations (Einasto 1965;Navarro et al. 2004;Wang et al. 2020).It can be written as where  0 and  s are the normalization and scale radius respectively, and  is the inverse of the Einasto index which characterizes the mass concentration.In this work, we only choose two parameter benchmarks for the Einasto profile as a representation of different Einasto index and local DM density.The first benchmark labelled as "EinN04" is originated from the simulation (Navarro et al. 2004): { ein,halo (8.2 kpc) = 0.01  ⊙ pc −3 ,  s = 20 kpc,  = 0.17}.However, the second one labelled as "EinW22" is derived by the rotation curve and globular cluster kinematics from Gaia (Wang et al. 2022): { ein,halo (8.2 kpc) = 0.008  ⊙ pc −3 ,  s = 12 kpc,  = 0.32}.We adopt the circular-orbit approximation of the GS spike model (Blumenthal et al. 1986;Ullio et al. 2001), which only differs from the semi-analytical one by no more than a factor of two (Ullio et al. 2001).The spike profile results from the conservation of the angular momentum, the radial action and the phase space distribution during the adiabatic growth of the BH (Young 1980).In the circularorbit approximation, the radial action is zero, the angular momentum  Figure 2. The logarithmic slope of the Einasto profile with and without the GS spike in solid and dashed lines respectively.The blue and orange lines correspond to the Einasto profiles with the parameter sets labeled as EinN04 (Navarro et al. 2004) and EinW22 (Wang et al. 2022), respectively.The inset shows the slopes of the GS spike profiles for the two models within 1 pc.The shaded region covers the orbit of S2.
conservation is given by and the conservation of distribution implies where  tot () and  dm () are the enclosed total mass and DM mass inside the radius  respectively, and the subscript  and  denote the initial and final state.The the mass distribution of DM before the growth of SMBH are Given the DM mass distribution after the growth  dm,  , the spike density profile can be calculated with where the inner radius 2 sch is for the Schwarzschild BH (Sadeghian et al. 2013;Ferrer et al. 2017).
We depict the Einasto profile and the GS spike profile by the blue solid and dashed lines in Fig. 1.Although the Einasto profile is smaller than NFW profile by two orders of magnitudes in the S2 orbit region, the spike profile after accretion is quite similar.The slopes of the GS spikes in the EinN04 and EinW22 benchmarks are presented by the blue and orange dashed lines in Fig. 2. The spike slopes barely change within the spike radius and are equal to 2.26 and 2.28 for the EinN04 and EinW22, respectively.This can be understood that the phase space distributions of the Einasto profiles are singular around the BH (Cardone et al. 2005;Baes 2022) so that the Einasto spikes are as steep as the gNFW spike with  = 0.
To release the GS assumptions, we choose the spike radius  sp as the free parameter.To obtain the enclosed mass profile given  sp , we use and where  acc ( sp ) is the mass of accreted interstellar medium, which is set to be  acc ( sp ) =  dm, (2 sp ).With Eq.( 6) and Eq.( 7), the enclosed mass profile  dm,  (;  sp ) and thereby the density profile  ein,sp (;  sp ) can be derived.If the spike radius  GS sp is used,  dm, (2 GS sp ) =  dm,  ( GS sp ) =  BH , therefore the enclosed mass profile  dm,  (;  GS sp ) can reduce to the GS spike mass distribution.

Dynamical model
To derive the stellar orbits, we solve the first-order post Newtonian approximation considering that the higher order effects are below the precision of current observation sensitivities (Do et al. 2019).The equation of motion with a spherically symmetric mass distribution can be written as (Rubilar & Eckart 2001) where  tot () =  BH +  dm () is the enclosed total mass which is related to the parameters of the spike.() = − ∫  ∞ d   tot ()/ 2 is the gravitational potential at a given radius.We define a coordinate system by setting the origin at the SMBH, letting the  and  axes point to the west and north, and making the  axis point from the GC to the solar system (Do et al. 2019;Yuan et al. 2022a).In this case, r() ≡ [ (),  (),  ())] and v() ≡ r() = [  (),   (),   ()].
An initial condition is also required to solve the equation of motion.We define the initial state for each star at the epoch  0 = 2000.0with the following six parameters: the inclination , the longitude of ascending node Ω, the positions ( 0 ,  0 ) and velocities ( 0 ,  0 ) in the orbital plane.The initial phase-space coordinates can be transformed from the six parameters through The DOP853 algorithm (Hairer et al. 1993) in SciPy is utilized to solve the state of the star at epoch .
The orbital planes of the S-stars are not perpendicular to the line of sight, so there is a time delay caused by the propagation of light through the orbit plane in  direction, i.e. the Rømer delay.This effect is detectable in the Keck observation of S2, which leads to a time delay of −0.5 days at pericentre and 7.5 days at apocentre (Do et al. 2019).The time delay is  obs =  em −  ( em )/, where  obs and  em represent the epochs of observation and emission respectively.The epoch of emission can be solved with iterations (Hees et al. 2014), however one iteration is adequate at the present (Do et al. 2019) The second term has an opposite sign compared to Do et al. (2019), since the  axis in their work is pointing from the Sun to the GC.The Shapiro time delay is ignored in our work since the correction is merely ≲ 5 min.
Once we have the star position at any epoch, the right ascension (R.A.), declination (Dec.) and radial velocity can be calculated with where  0 is the distance between the Sun and the central SMBH and  ref is the reference epoch which is 2000.0 for Keck (Do et al. 2019) and 2009.0 for VLT (Gillessen et al. 2017).( BH ,  BH ) and ( ,BH ,  ,BH ) are the offset and linear drift of the central mass in celestial coordinate. = |v| is the norm of the velocity, and  0 is a constant velocity offset in the radial velocity measurements.

Data and statistical method
We use the latest publicly available astrometric and spectroscopic measurements of the S-stars from the VLT and Keck telescopes compiled in Do et al. (2019) and Gillessen et al. (2017).Do et al. (2019) presents the measurements of S2 from Keck/NIRSPEC, Keck/NIRC, Keck/NIRC2, Keck/OSIRIS, Gemini/NIFS, Subaru/IRCS and VLT/SINFONI between 1995 and 2019.Gillessen et al. (2017) provides the latest publicly available measurements of 17 S-stars observed by NTT/SHARP, VLT/SINFONI, VLT/SPIFFI and VLT/NACO from 1992 to 2017.Since the radial velocities compiled in Do et al. (2019) contain the VLT/SINFONI measurements in Gillessen et al. (2017), we remove the overlaps in the combined analysis.Different offsets and drifts of the central object are adopted as free parameters for two observatories to explain their discrepancies in the results (Gillessen et al. 2009;Abuter et al. 2021b), which enables us to combine the data sets.For the S2, we have 190 astrometric points and 115 radial velocities with a full coverage of the orbit.Although there are updated measurements of S2, S29, S38 and S55 from the VLT/GRAVITY interferometric and VLT/SINFONI spectroscopic observations (Abuter et al. 2018a(Abuter et al. , 2020(Abuter et al. , 2022)), these data are not publicly available yet.
The total likelihood function is constructed as follows: where L astro,Keck and L astro,VLT are the likelihood functions of the astrometric data from Keck and VLT respectively, L rv is the likelihood function for the radial velocity.The total likelihood contains the following free parameters: (i) Two parameters concerning the SMBH mass  BH and the distance to the GC  0 .
(ii) The offset ( BH ,  BH ) and the drift ( ,BH ,  ,BH ) of the central object in the celestial coordinate for each telescope.There are eight free parameters.
(iii) The velocity offset  0 and the additional systematic offset  offset between the Keck/NICR2 and other instruments (Chu et al. 2018;Do et al. 2019).
(v) One parameter on the DM spike: either the spike radius  sp or the density slope of the spike  sp .
When  stars are involved, there are 15 + 6 free parameters in the likelihood function.
For the astrometric data from the Keck, the correlation caused by the faint source confusion (Plewa & Sari 2018) is considered as recommended in Do et al. (2019).Therefore the likelihood function of the astrometric data from Keck reads: where Δ ≡ {  −   (  )} and Δ ≡ {  −   (  )} are the vector differences between the observed and predicted astrometric data from the Keck.
considering different noise models only change the parameters within the statistical uncertainty for the data set (Abuter et al. 2021b).
The likelihood function for the radial velocities is where   , and    (  ) are observed and predicted radial velocity, and    , is the uncertainty.The Bayesian inference method is chosen in this work.Two different approaches are adopted to analyze the spike profile.In the first approach, we fix the spike slope  sp to the value in the GS model and constrain the spike radius  sp .We assume a flat prior for  sp over a range greater than 2 sch (Lacroix 2018) in this case.In the second approach, we fix the spike radius  sp to the GS predicted value and constrain the spike slope  sp .A prior distribution for  sp is assumed based on the flat prior of the enclosed DM halo mass within the S2 apocentre  apo , namely where  apo ≡  apo / sp and  sp < 3. The priors of the remaining parameters are set to be flat in linear space.We adopt the Markov chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013) to sample the posterior probability.To ensure the samples well converged, we require the Gelman-Rubin diagnostic of  − 1 ≲ 0.05 and the samples longer than ∼ 30 times the integrated autocorrelation per chain (Hogg & Foreman-Mackey 2018;Cowles & Carlin 1996).

Constraint on the generalized NFW spike with S2
The gNFW spike profile contains two parameters, the spike radius  sp and the spike slope  sp , which can be altered due to the massive seed BH or the dynamical heat process.Since the S2 orbit only covers a small range of the spike (the dark grey band in Fig. 1), it is difficult to determine the two parameters at the same time.Hence, in this subsection, we treat the upper limits of  sp and  sp separately, and keep the other parameter fixed to the value given in the GS model.

The spike radius 𝑅 sp
We analyze the combined orbital data of S2 with the gNFW halo model including a spike.In this part, we constrain the spike radius  sp and fix the spike slope to the GS predicted value.In Fig. 3, we plot the measurements of S2 overlaid with the maximum-a-posteriori (MAP) spike model with  = 1.The left panel is the celestial coordinates from the SMBH corrected for the offsets and drifts of the two telescopes.The orange hollow and blue solid points are the astrometric measurements from the VLT and Keck telescopes respectively, while the black line represents the MAP orbit of S2 which moves clockwise.The top, middle and bottom panels on the right are the R.A., Dec. offsets and radial velocity as a function of time.The MAP model can well fit the measurements.
Following Lacroix (2018), we take into account different gNFW density slopes  and analyze the constraint on the spike radius  sp for each case.In the left panel of Fig. 4, we show the 68% credible interval and 95% upper limit of  sp with the colour band and solid line.The spike spatial extensions larger than 71.9 pc, 15.7 pc and 0.56 pc for the gNFW slope  of 0.5, 1.0 and 1.5 are excluded at 95% level.The  sp constraint for a steeper initial slope  tends to be more stringent since it causes more DM mass within the range of S2 orbit.We also draw the 95% upper bounds from Lacroix (2018) with a pink dot-dashed line in the figure.Our constraints are much stronger thanks to accurate Keck measurements of S2 at the closest approach.The predicted spike radius in the GS model (Eq.( 4)) is shown in the black dashed line and the ones with  ≥ 0.92 are excluded at 95% level.Even when we consider the systematic uncertainty of the local DM density  halo,⊙ = 0.008 − 0.013  ⊙ pc −3 (de Salas & Widmark 2021), the predicted spike radius for the initial gNFW with  = 1 ranges from 15.8 pc to 20.5 pc and is still inconsistent with data.
We also convert the  sp to the enclosed DM spike mass within the apocentre of S2.As shown in the right panel of Fig. 4, the constraint only weakly depends on , confirming that the combined orbital data of S2 is not sensitive to the slope of the distribution at the present (Do et al. 2019).The 95% credible upper limit of enclosed mass is found to be ∼ 6000  ⊙ , when the combined S2 data are adopted.The extended mass constraint is weaker than that of the GRAVITY (Abuter et al. 2020), since the latter contains the private interferometric measurements which are more accurate than the adaptive optics.
The constrained spike models are consistent with other types of observations.Such a spike increases the rotation velocity of gas by ≲ 2% within the Galactocentric radius ≲ 5 pc, which is much smaller than the statistical uncertainties of the rotation curve (15% − 30%) (Sofue 2013).The spike contributes to an extended mass of ≲ 20  ⊙ within 2 × 10 −6 pc, which is not in conflict with the measurements of three hot spots around the SMBH considering the uncertainty of the inferred enclosed mass of ∼ 10 6  ⊙ (Abuter et al. 2018b).The median precession angle of S2 is 10 ′ − 11 ′ per orbit when such a spike exists (Fig. A1), consistent with the GRAVITY measurements within ∼ 1 uncertainty (Abuter et al. 2020).

The spike slope 𝛾 sp
We analyze the spike slope  sp using the measurements of S2.The spike radius is fixed to the predicted value in the GS model.The same as the previous part, we test the gNFW models with different halo density slope .In Fig. 5, we present the 68% credible interval and 95% upper limit of  sp with the colour band and solid line.The spike slope  sp steeper than 2.45, 2.32 and 2.14 for the  of 0.5, 1.0 and 1.5 are disfavoured by the S2 orbit at 95% probability.We draw the spike slopes predicted by the GS spike model (Eq.( 3)) in the black dashed line.The GS spike with  ≥ 0.92 is excluded at 95% level.
The DM spike slope  sp could be flattened to 3/2 when the DM particles are efficiently heated by the dynamical processes (Merritt et al. 2002;Gnedin & Primack 2004;Shapiro & Heggie 2022).The DM spike profile may also resemble the star spike profile in the GC whose slope  sp is 1.1 − 1.6 (Schödel et al. 2020;Abuter et al. 2020), considering the collisionless nature of DM particles.However, it is difficult to constrain such possibilities with the current S2 measurements at the present.

Constraint on the Einasto spike with S2
The Einasto profile is among the most popular models for the DM density profile, however, the corresponding spike profile has not been discussed and constrained in the literature yet.Here we set upper limits on the spike radius in the two parameter benchmarks of the Einasto spike profile given in Sec.2.2.Fig. 6 presents the marginal posterior probability of spike radius  sp for two Einasto parameter sets.The dark and light shaded regions represent the 68% and 95% credible upper limits.For the benchmark EinN04 and EinW22, the 95% upper limits of the spike radius are 21.5 pc and 61.4 pc respectively.Comparing to the expected spike radii of 31.2 pc and 60.4 pc in the GS model shown in the black dashed lines, the steep EinN04 spike is disfavored at 95% level and the flat EinW22 spike marginally survives from the S2 constraint.

Joint analysis with multiple S-stars
Gillessen et al. ( 2017) has reported the VLT observations of an additional 16 S-stars.The combination of multiple S-stars can span a broader range of density profile, as illustrated by the shaded band in Fig. 1, thereby potentially improving the constraint on the spike density profile.This subsection is to discuss the outcomes obtained when multiple S-stars are combined.
The VLT data of S1, S9, and S13 are incorporated into the previous S2 dataset, as suggested by Gillessen et al. (2017) due to their increased sensitivity in constraining  0 and  BH .It can potentially distinguish the extended mass from SMBH, and provide a stringent constraint on the DM spike.Including 637 astrometric points and 228 radial velocities, the constraints on the gNFW spike are re-evaluated.The  2 values of the three stars are combined with Eq.( 18) and Eq.( 19), resulting in a total likelihood function containing 39 free parameters.
In Fig. 7, we show the measurements of the four S-stars (points) overlaid with MAP models (solid lines).The left panel is the celestial coordinates corrected for the offsets and drifts.The right panel is the radial velocities with respect to the epoch for the stars.Indeed, we can see that the MAP orbital models can match the measurements reasonably well.
In Fig. 8, we present the upper limits on the spike radius  sp with the four S-stars.The constraints from the combined analysis are slightly stronger than the case only including S2 data.The expected spike radius in the GS model is shown in the black dashed line.We show that the GS model with  ≥ 0.90 is excluded by the observations of the four stars at 95% probability, by comparing with  ≥ 0.92 in the case of S2 only.
In Abuter et al. (2022), the GRAVITY data of S2, S29, S38, and S55 are combined to further improve the constraints of extended mass.We take the publicly available VLT data of S38 and S55 (Gillessen et al. 2017) and incorporate them into the previous  S2 data set.In total, we have 348 astrometric points and 122 radial velocities.Similar joint analyses are performed and the constraints on  sp are presented in Fig. 8.The constraints are slightly stronger than the results using four S-stars, due to the small uncertainties of their semi-major axes and orbital periods.The GS model with  ≥ 0.83 is excluded at 95% probability.

Interpretations of GRAVITY present and future S2 sensitivity of the halo extended mass on DM spike
The GRAVITY coherently combines the light of the Very Large Telescope Interferometer and presents ∼ 10 as-level photometric observations of S-stars (Abuter et al. 2017).In Abuter et al. (2022), the S2 data collected from 1992.2 to 2021.6 were adopted to constrain the extended mass.They have found that no more than 2400 ⊙ or 7500 ⊙ extended mass could exist within the apocentre of S2 at 1 or 3 level respectively.In addition, Heißel et al. (2022) have estimated a prospect that the GRAVITY could gather a full-orbit data of S2 with 50 as and 10 km s −1 precision by 2033.Their result shows that the future S2 sensitivities from GRAVITY could restrict the dark extended mass less than 1000 ⊙ at 1 level.Hence, in this section we simply project the GRAVITY extended mass limits (both present combined limits and the future sensitivities) on  sp and  sp .
We learn from the left panel of Fig. 4 that the enclosed mass within the apocentre of S2 only weakly depends on the density slope.We can use the extended mass constraints derived from the Plummer or Bahcall-Wolf cusp profiles in Abuter et al. (2022) and Heißel et al. (2022) as approximated upper limits of the extended mass of the DM spike, even though the slopes of the two profiles are different from the DM spike.In this way, we directly convert the extended mass upper limits to the constraints of DM spike parameters.
The left panel of Fig. 9 shows the upper limits on the spike radius  sp for the gNFW spike.The purple solid and dashed lines corre- Figure 9.The upper limits (ULs) on the spike radius  sp (left) and slope  sp (right) of the gNFW spike profile converted from the GRAVITY upper limits of the extended mass within the apocentre of S2.The solid and dot-dashed purple lines come from the 1 and 3 constraints up to 2021 in (Abuter et al. 2022), while the brown solid line is from the expected 1 constraint from a full orbit of the GRAVITY data (Heißel et al. 2022).The expected spike radius in GS spike model is shown with the black dashed line.
spond to the 1 and 3 upper limts on  sp by using the GRAVITY data collected before 2021.6.For the DM halo with  = 1, the 1 and 3 upper limits are 8.0 pc and 19 pc, respectively.Our results are in good agreement with the ones obtained by fitting the GRAVITY orbital data ( sp ≲ 10 pc) (Abuter et al. 2020).The black dashed line in the figure shows the corresponding  sp in the GS spike model [Eq.( 4)].The gNFW spike with the initial density slope  steeper than 0.64/1.00 is rejected by the data of the current GRAVITY S2 measurements at 1/3 level.The future sensitivities of  sp are presented by the brown solid lines in the case of using a full-orbit GRAVITY data of S2.The 1 upper limit on  sp for the gNFW with  = 1 is ∼ 4 pc in this case, and the GS spike with  ≳ 0.4 can be probed.
The right panel of Fig. 9 illustrates the limits of the spike slope  sp for the gNFW spike.Given the GRAVITY measurements collected before 2021.6, the 1 (3) upper limits on  sp for the gNFW with  = 1 are 2.21 (2.34).Once a full-orbit GRAVITY data of S2 were available, the prospective 1 limit could be improved, and the GS spike with  sp > 2.10 can be completely probed.With the GRAVITY S2 orbit only, it is far to detect the DM spike slope  sp ∼ 1.5, corresponding to the DM particles efficiently heated by the dynamical processes.
We also estimate the GRAVITY constraining power for the Einasto spike.The 1 (3) upper limits on  sp for the EinN04 and EinW22 benchmarks are 15 pc (30 pc) and 45 pc (83 pc), respectively.The GS model for the EinN04 benchmark ( GS sp = 31.2pc) is excluded at 3 level.If adopting the full-orbit GRAVITY measurements, the prospective 1 upper limits could be strengthened to 9 pc and 29 pc for EinN04 and EinW22, respectively.The GS Einasto spike model can be entirely probed by the future full-orbit GRAVITY data of S2.

The effect of dark matter annihilation
The DM spike distribution is widely used in the DM indirect detection, DM particles could annihilate with each other and produce detectable signals.In previous sections, we have discussed the cases of lacking DM annihilation.Since the DM annihilation can reduce the spike density, as long as the annihilation cross section ⟨⟩ is large enough, the DM spike may also survive from the S2 orbit con-straints.In this subsection, we present such a lower limit on ⟨⟩ for the NFW spike.
Considering that the initial velocity distribution of DM particles is isotropic, the particles with apicentres outside  sat can also contribute to the density inside  sat , leading to a weak cusp with a slope of 0.5 for the -wave annihilation (Vasiliev 2007;Shapiro & Shelton 2016).The density profile below is adopted (Vasiliev 2007) where  sat ≡  DM /(⟨⟩  BH ) is the saturation density,  sat is the saturation radius defined by  GS sp ( sat ) =  sat ,  DM is the DM particle mass, and  BH = 10 Gyr is the age of the SMBH in the GC (Lacroix et al. 2014).The NFW spike density distribution and the spike radius in the GS model are  GS sp () and  GS sp , respectively.We choose the saturation density  sat as the free parameter, set a flat prior for the enclosed mass within the apocentre of S2, and rerun the MCMC sampler using the combined Keck and VLT data.By marginalizing posterior probability, the 95% credible upper limit of  sat is ∼ 2.2 × 10 9  ⊙ pc −3 .The corresponding saturation radius is  sat ∼ 4.8 mpc.Therefore, the surviving NFW spike infers the 95% lower limit of the annihilation cross section: ⟨⟩ ≳ 7.7 × 10 −27 cm 3 s −1 ×  DM 100 GeV 10 Gyr Interestingly, such a request is indeed satisfied by a good fraction of the parameter space of the inert two Higgs doublet model for the W-boson mass excess (Fan et al. 2022).

SUMMARY
The density of DM can be significantly increased due to the adiabatic growth of BH, creating a structure so-called DM spike.In this work, we constrain the parameters of the DM spike profile in the GC by using the precise measurements of the S-stars from the Keck and VLT observatories.Table B1.The 95% credible upper limits on the spike radius  sp or spike slope  sp for different initial NFW density slope .The values in the second and fifth columns are the values predicted by Gondolo & Silk (1999).The third and sixth columns show the constraints using the VLT and Keck data of S2.
The fourth column presents the constraints using the combined observations of 3 S-stars.See Sec. 4 and Sec.5.1 for detail.

Figure 1 .
Figure 1.The possible DM density profile in the Milky Way.The red and blue lines are for the initial NFW and Einasto profile, respectively.The solid and dashed lines represent the profiles with and without the GS spikes.The radial extension radii are also drawn in dots.The ranges covered by S2-only and the four stars (S2, S1, S9 and S13) are indicated, respectively, by the dark shaded region and the entire shaded region.

Figure 3 .Figure 4 .Figure 5 .
Figure3.The data and the maximum-a-posteriori (MAP) model of S2 orbit around the SMBH with the DM spike of  = 1 included ( sp is free).In the left panel, astrometric measurements corrected for the offsets of the reference frames observed by the VLT (Keck) are marked with orange hollow (blue solid) points.The black solid line shows the MAP orbit with the S2 moving clockwise.The Sgr A* is located at the origin marked with a black star.The right ascensions, declinations and radial velocities as a function of the epoch are shown in the upper, middle and lower right panels.The measurements are marked with points, while the optimal models are drawn with the dashed lines.Since there is a velocity offset for the Keck/NIRC2 data comparing to the rest data, we plot them separately in green square points.

Figure 6 .
Figure 6.The marginal posterior distribution of the spike radius  sp of the Einasto spike.The left and right panels correspond to the Einasto spike profiles labeled as EinN04 and EinW22.The dark and light shaded regions represent the 68% and 95% credible upper limits.The dashed line shows the spike radius which results from the adiabatic growth of the SMBH  BH as given in Sec. 2. The spike profile initiated from the profile EinN04 is disfavored by 95% probability.

Figure 7 .
Figure7.The measurements and the MAP models of S2, S1, S9 and S13 around the SMBH with the DM spike of  = 1 included ( sp is free).The left panel shows the astrometric points corrected for the offsets of reference frames, while the right panel shows the radial velocities for the epoch of observation.The hollow and solid points represents the data from the VLT and Keck respectively.The solid lines illustrate the MAP orbits and radial velocities.S2 moves clockwise, while S1, S9 and S13 move anti-clockwise.

Figure A1 .
Figure A1.The S2 precession angle per orbit Δ when DM spike with initial NFW slope  exists around the SMBH.The blue solid line and blue band correspond to the median and 68% uncertainties of the precession angle respectively.The fitted value and the uncertainty measured by the GRAVITY (Abuter et al. 2020) are shown with the black solid line and black band, respectively.The orange dashed line illustrates the Schwarzschild precession angle.
(Gillessen et al. 2017) are [  ]   ≡          and [  ]   ≡          , where   and   are the uncertainties of the astrometric data, and  is the correlation matrix defined as []  = (1 − )   +  exp −   /Λ ,where    is the angular distance between the two Keck data points  and .We also set free the correlation length scale Λ and the mixing parameter  in the model.On the other hand, we simply use the  2 form of the likelihood function for the VLT astrometric data(Gillessen et al. 2017)−2 ln L astro,VLT = ∑︁ −   (  )   2 sch ,  sat / √︁ / sat 2 sch <  ≤  sat ,