Dynamical constraints on a dark matter spike at the Galactic Centre from stellar orbits

In this work I use astrometric and spectroscopic data on the S2 star at the Galactic Centre (GC) up to 2016 to derive specific constraints on the size of a dark matter (DM) spike around the central supermassive black hole Sgr A*. These limits are the best direct constraints on a DM spike at the GC for non-annihilating DM and exclude a spike with radius greater than a few tens of parsecs for cuspy outer halos and a few hundred parsecs for cored outer halos.


Introduction
Dark matter profiles in the central regions of galaxies are poorly constrained at present and are the objects of intense debate. While observations seem to favour flat (cored) profiles, numerical simulations favour steeper profiles (cusps), leading to the cusp/core controversy (e.g. de Blok 2010 for a review). At subparsec scales, the dark matter (DM) distribution is even less constrained and can be significantly affected by the central supermassive black hole (SMBH). In particular, if the SMBH grows adiabatically, i.e. on a much longer timescale than the dynamical timescale, the DM density is expected to be significantly enhanced (by up to 10 orders of magnitude at the very centre) in a region corresponding to the sphere of influence of the black hole (BH), typically at parsec scales for the Milky Way. This leads to a very sharp morphological feature referred to as a DM spike, corresponding to a DM profile going as r −γ sp , with γ sp typically between 2.25 and 2.5, depending on the slope of the initial DM halo (Gondolo & Silk 1999). DM spikes are of particular interest in the context of indirect DM searches since they lead to very strong signatures of DM annihilation and allow us to probe weakly annihilating DM particles (Gondolo & Silk 1999;Regis & Ullio 2008;Lacroix et al. 2014Lacroix et al. , 2015Lacroix et al. , 2017Fields et al. 2014;Shapiro & Shelton 2016).
There is, however, considerable uncertainty on the formation and survival of DM spikes. In particular, the assumption of adiabaticity may not be verified in general. For instance, dynamical processes such as mergers can lead to weaker cusps (Merritt et al. 2002). However, binary scouring only occurs above parsec scales, while we are interested in the DM profile much closer in when studying the orbits of S stars, as discussed in the following. Moreover, the Milky Way (MW) is unlikely to have suffered such mergers in its recent past, as evidenced by the quiet history of the thick disk since the only major merger which occurred about 12 Gyr ago and is likely to have led to the formation of the bulge and the SMBH (Wyse 2001). A weaker cusp is also formed if the BH does not grow exactly at the centre of the DM halo (within ∼ 50 pc) (Nakano & Makino 1999;Ullio et al. 2001) or if the BH growth cannot be considered adiabatic Ullio et al. (2001), but the actual impact of these effects on the MW is unclear. Moreover, dynamical heating in the central stellar core would also soften a spike (Gnedin & Primack 2004). Another concern is that the non-observation of a stellar spike (recent results point to a softer stellar cusp than previously thought, with slope ∼ 1.15; Schödel et al. 2017) would rule out the existence of a DM spike. However, if the BH grows (for example by gas accretion) mostly before the nuclear star cluster forms in the spike region, then the DM and stellar profiles are decoupled. Additionally, the nuclear star cluster in the most accepted view is formed by merging globular clusters, as in Antonini et al. (2015). This leads to different profiles for the DM and stellar distributions. Therefore, stars and DM essentially decouple, and the absence of a stellar spike in observations does not preclude the existence of a DM spike. On the other hand, additional dynamical processes can have the opposite effect of regenerating a spike, for example enhanced accretion of DM to counteract the depopulation of chaotic orbits in triaxial halos (Merritt & Poon 2004) or gravo-thermal collapse for self-interacting DM (Ostriker 2000).
As a result the unclear status of the inner DM profile of galaxies as discussed above calls for direct probes. In particular, there is still no definitive evidence either in favour of or against such a high concentration of DM either in the MW or in any other galaxy. This is due in particular to the small size of the regions involved. Probing such regions requires high angular resolution and astrometric precision to characterize the gravitational potential. However, the inner region of the MW offers a unique window on the DM distribution at the Galactic Centre (GC), thanks to the monitoring of the orbits of the S stars within ∼ 1 arcsec of the central BH. In particular, since it is the closest star to the BH observed so far, the S2 star has been extensively studied through monitoring campaigns based on observations conducted with the Very Large Telescope (VLT) (Schödel et al. 2002;Gillessen et al. 2009aGillessen et al. ,b, 2017GRAVITY Collaboration 2018) and the Keck observatory (Ghez et al. , 2008 Article number, page 1 of 8 arXiv:1801.01308v2 [astro-ph.GA] 4 Sep 2018 Boehle et al. 2016). 1 These series of observations have led to the reconstruction of the orbit of the star over roughly one and a half periods. In addition to tight constraints on the mass of the central SMBH, M BH , and its distance from Earth, R 0 , these two groups have shown that only a small fraction (typically 1-2%) of the mass of the SMBH can be in the form of an extended distribution (Ghez et al. 2008;Gillessen et al. 2009a,b;Boehle et al. 2016;Gillessen et al. 2017). Other constraints have been obtained on an extended component by studying the corresponding reconstructed mass profile (Hall & Gondolo 2006) or the pericentre shift of S2 (Zakharov et al. 2007;Iorio 2013).
Here I go a step further and I use astrometric and spectroscopic measurements of the orbit of S2 up to 2016 to set specific constraints on the DM distribution in the inner Galaxy. I present the first direct dynamical constraints from stellar orbits on the size of a DM spike, inside a DM halo constrained by larger scale kinematic data at kpc scales, for example from maser observations. This is especially interesting for non-annihilating or very weakly annihilating DM which is not expected to have significant observational signatures other than gravitational.
In Sect. 2 I describe the model along with the orbit-fitting procedure, before presenting my results in Sect. 3. Finally, I conclude in Sect. 4.

Model and orbit-fitting procedure
2.1. Calibration: the point-mass case I rely on textbook results of standard mechanics in a central potential (e.g. Bate et al. 1971). I first recall the parameters of the problem in the BH-only case, which has an analytic solution, before moving on to the more general case of an extended mass distribution. The BH-only case serves as calibration for the orbitfitting procedure.
The orbit-fitting procedure consists in reconstructing the time evolution of the position and velocity of the star on its orbit to determine the properties of the gravitational potential by fitting the parameters of the model to the data. In the case of one star orbiting a central point mass, the 13 parameters of the problem are the mass of the central object, here denoted M BH , and its six phase-space coordinates, namely its distance R 0 , its position on the sky (α BH , δ BH ), and velocity (v α,BH , v δ,BH , v r,BH ), as well as the six phase-space coordinates of the star. However, the orbit of the star is more readily characterized analytically in terms of the six standard orbital elements: the semi-major axis a of the orbit, the eccentricity e, the time of pericentre passage t P , and three angles, namely the inclination I of the orbital plane with respect to the plane of the sky, the longitude of the ascending node Ω, and the angle ω between the directions of the ascending node and the pericentre.
Although the motion of Sgr A* with respect to the local standard of rest (LSR), defined as the circular velocity at the radius of the Sun, is expected to be very small (Reid & Brunthaler 2004;Plewa et al. 2015), its position (α BH , δ BH ) on the plane of the sky at a reference time t ref and its velocity (v α,BH , v δ,BH , v r,BH ) relative to the LSR are unknown a priori and can be constrained through the orbit-fitting procedure. In practice, the motion of the BH is accounted for through a linear term in the angular position of the star as a function of time. The reference time is taken to be 2009 yr (Gillessen et al. 2017) for the data set up to 2016, and 2005.4 yr for the data set up to 2009 (Gillessen et al. 2009b).
1 Data from before 2002 were produced with the New Technology Telescope (NTT).
In this work, I used the data from the NTT/VLT and Keck observatories compiled in Boehle et al. (2016); Gillessen et al. (2017). In Gillessen et al. (2009b) the authors presented a robust method to consistently combine the two independent data sets for which the astrometric data feature a clear offset due to slight differences in the definition of the coordinate systems. More specifically, to account for the discrepancy between the two data sets, they introduced an offset in angular position (∆α, ∆δ) and velocity (∆v α , ∆v δ ) on the plane of the sky to shift the Keck data back onto the VLT data. This was done by fitting the model with these 4 parameters in addition to the 13 parameters described before. In practice, this is achieved by shifting the observed right ascensions and declinations of the star measured with the Keck observatory by the quantities ∆α + ∆v α (t − t ref ) and ∆δ + ∆v δ (t − t ref ), respectively. I repeated this procedure here for the combined data set.
Throughout this work, I derived the posterior probability density function of model parameters using PyMultiNest (Buchner et al. 2014), which relies on the MultiNest code (Feroz et al. 2009) based on the multimodal nested sampling Monte Carlo technique (Feroz & Hobson 2008). Multimodal nested sampling is particularly suitable for studying high-dimensional parameter spaces with possible degeneracies between parameters. The likelihood combines the data on right ascension, declination, and radial velocity of S2. I used uniform priors for all parameters except the position and velocity of the BH, for which I took Gaussian priors based on the results from Plewa et al. . I also recovered the best-fit model for the combined VLT+Keck data set, using the prescriptions of Gillessen et al. (2009b) for the priors on ∆α, ∆δ, ∆v α , and ∆v δ . This served as a consistency check of the analysis chain, which was then applied to the study of the impact of a DM spike on the orbit of S2.

Extended mass
For the extended DM mass distribution, I consider two scenarios: the general case of a non-annihilating cold dark matter (CDM) candidate, and the more specific case of self-annihilating DM, applicable to candidates like weakly interacting massive particles (WIMPs). For non-annihilating DM the spike goes way inside the orbit of S2, down to the close vicinity of the SMBH (Sadeghian et al. 2013): where R sp is the radial extension of the spike, and the halo profile is assumed to be given by a generalized Navarro-Frenk-White (NFW) profile characterized by a slope index γ, where r s is the scale radius, and the scale density ρ s is related to the local density ρ via (2) More specifically, the idea is to consider the various DM halos corresponding to the dynamically constrained Milky Way mass models from the analysis of McMillan (2017), and determine the maximum size of a DM spike inside that halo that does not cause a significant departure from the best-fitting BH-only orbit. The associated values of the local density, scale radius, and R 0 from the analysis of McMillan (2017) are summarized in Table A.1, in Appendix A. 2 For self-annihilating DM, the inner region of the DM spike is depleted since the DM density is so high that DM particles annihilate more efficiently. This results in a plateau of density ρ sat = m DM /( σv t BH ), where m DM is the mass of the DM candidate, σv is the velocity-averaged annihilation cross section, and t BH is the age of the central SMBH, which I take conservatively to be ∼ 10 10 yr. The saturation plateau extends to a radius The cases of non-annihilating and self-annihilating DM are both illustrated in Fig. 1. Shown are the density profiles for regular NFW-like halos and halos with a spike in the central region for non-annihilating DM and a self-annihilating 1 TeV DM candidate with three values of σv (see figure for details). The corresponding mass profiles are shown in the right panel of Fig. 1 with the same line styles. The profiles are illustrated with γ = 1, which gives a spike slope γ sp = 7/3, and a spike radius R sp ∼ 100 pc which corresponds to the 99.7 % upper limit I obtain from deviations of the BH-only orbit, as discussed in Sec. 3. For the self-annihilating case, the values of the cross section are chosen to illustrate the point at which the annihilation plateau becomes as big as the characteristic size of the or-bit, given by the semi-major axis constrained to be of the order of 5 mpc by the orbit-fitting procedure. For m DM ∼ 1 TeV and σv ∼ 10 −26 cm 3 s −1 , the mass enclosed inside the orbit is significantly reduced with respect to the case of non-annihilating or very weakly annihilating DM, down to values much smaller than a few percent of the BH mass, making deviations from the BH-only orbit undetectable. 3 In the absence of a spike, the DM halo has a negligible impact on the orbit of S2, as illustrated in the right panel of Fig. 1, due to a much smaller mass enclosed in the orbit. Moreover, for completeness I have also considered the effect of a realistic stellar profile ρ star ∝ r −γ star , with γ star ∼ 1.15 and ρ star (1 pc) = 1.5 × 10 5 M pc −3 (Schödel et al. 2017). However, the corresponding mass enclosed in the orbit is about three orders of magnitude below the critical mass needed to have an impact on the orbit. The same conclusion applies to the stellar bulge profile from McMillan (2017).
In the general case of an extended mass distribution around the central point mass, the orbit model is no longer analytic and one must rely on numerical tools to solve the equations of motion. First the polar radius r(t) is determined with Newton's second law in the Galilean frame of the BH, where L ≡ r 2θ is the angular momentum modulus, G is the gravitational constant, Φ ext is the potential created by the extended mass, and Φ S accounts for the effect of Schwarzschild precession induced by the BH. 4 The initial conditions r 0 ≡ r(t 0 ) and 3 These considerations can be extended to other values of m DM via Eq. (3). 4 Although the current data on S2 are not yet sensitive to relativistic effects (Gillessen et al. 2017), I include this post-Newtonian precession effect for completeness since it is partly degenerate with the precession caused by the extended mass and as such can mildly affect the limits set on the DM profile. Other relativistic effects such as gravitational redshift essentially affect radial velocities, which are much less tightly constrained by observations than the position of the star, as discussed in Gillessen et al. (2009a).
A&A proofs: manuscript no. S2_final_arXiv  (2017), and linearly interpolated to get smooth curves. The prediction from Gondolo & Silk (1999) is also shown as a benchmark model (cyan solid).
r 0 ≡ṙ(t 0 ) need to be specified, where t 0 is chosen as the first epoch in the data, namely t 0 = 1992.224 yr. For given values of r 0 andṙ 0 , I use the odeint Python routine to solve for r(t). Once r(t) is known, θ(t) is obtained via where θ 0 ≡ θ(t 0 ) and L = r 2 0θ 0 , withθ 0 ≡θ(t 0 ). Orbital elements no longer characterize the orbit but only an osculating orbit in the general case of an extended mass, so they cannot be used to parametrize the problem. The free parameters for the star are now the initial conditions r 0 , θ 0 ,ṙ 0 ,θ 0 , as well as I and Ω, which still characterize the plane of the orbit of the star. The parameters of the BH do not change. 5

Results
In this section, I present constraints on the size of a DM spike as a function of the slope of the DM halo obtained with the multimodal nested sampling analysis implemented in PyMultinest. It should be noted that the nested sampling procedure does find a non-zero best-fit value for the spike radius R sp for all values of γ. However, the very mild increase in Bayesian evidence when adding a DM spike, which remains smaller than ∆ ln Z ≈ 3, is insufficient to claim any preference for the BH+spike model (Kass & Raftery 1995). As a result, the data are consistent with the BH-only model. Nevertheless, it is still possible to exclude large values of the spike radius that would lead to a large DM mass inside the orbit, and thus to large deviations of the orbit of S2.
The resulting 95% and 99.7% confidence contours in the γ-R sp plane are shown in Fig. 2  In order to keep the same coordinate system as in the BH-only case, I fixed the value of ω, which is no longer a free parameter since the pericentre is not defined in general, to the best-fit value in the BH-only case. 0, 0.25, 0.5, 0.75, 1, 1.25, 1.5. I interpolated the results to obtain a smooth limit. 6 Using the combined VLT+Keck data set up to 2016, I exclude at the 99.7% confidence level a DM spike with a spatial extension larger than 90 pc for an outer halo with γ = 1, and larger than 6 pc for an outer halo with γ = 1.5. For the combined 2009 and the VLT-only 2016 data sets, the limits are about a factor of 2 weaker. The analysis of McMillan (2017) seems to favour cuspy halos (γ ∼ 1), whereas a recent study on the dynamics of the Galactic bar favours a DM halo with slope γ < 0.6 (Portail et al. 2017). For such cored halos, our constraints are weaker, with a maximum spike radius of a few hundred pc. Nevertheless, these limits are the first direct constraints on a DM spike at the GC, valid for non-annihilating DM. This is especially interesting because it is applicable to any CDM candidate with no significant annihilation cross section.
It should be noted that the limit already improves by about a factor of 2 when going from the combined 2009 data set to the combined 2016 one. Significant additional improvements can thus be expected on the constraints with the more recent S2 data.
Constraints on an extended mass component can also be expressed more generally in terms of the total extended mass M ext inside the characteristic size of the orbit. At a 99.7% confidence level, for the combined data up to 2009, I find M ext 2×10 5 M , corresponding to ∼ 4%M BH , while for the complete data set up to 2016, I obtain M ext 4-5 × 10 4 M , corresponding to ∼ 1%M BH . These results are consistent with the upper limits on a general extended mass component from Gillessen et al. (2009bGillessen et al. ( , 2017. These limits are illustrated in the right panel of Fig. 1 by the horizontal dotted and solid black lines. For self-annihilating DM, the constraints derived here are valid for σv < 10 −30 cm 3 s −1 for a benchmark particle mass m DM = 1 TeV. As shown in the right panel of Fig. 1, for the same candidate mass, the mass enclosed in a sphere of radius of the characteristic size of the orbit of S2-typically 5 mpcis decreased by about a factor of 2 for σv ∼ 10 −27 cm 3 s −1 and by about a factor of 20 for σv ∼ 10 −26 cm 3 s −1 with respect to the mass in the absence of annihilation. As a result, the upper limits on R sp are weakened by about a factor of 10 for σv ∼ 10 −27 cm 3 s −1 , while no constraints can be set for σv ∼ 10 −26 cm 3 s −1 . More generally, when the radius R sat of the saturation plateau due to self-annihilations (see Eq. (3)) reaches the size of the orbit of S2, the DM mass enclosed inside the orbit becomes too small to induce significant deviations from the BH-only orbit.
For illustration, the prescription from Gondolo & Silk (1999) is also shown in Fig. 2 (cyan solid line). This corresponds to a spike radius defined by with α γ ≈ 0.293γ 4/9 , which only differs from ∼ 0.1 for γ 1. It should be noted that these predictions are indicative and can be significantly affected by the various dynamical processes discussed in Sec. 1. For very cuspy halos (γ ≈ 1.5), the combined 2016 data already exclude the prediction from Gondolo & Silk (1999) at a 95% confidence level.

Conclusion and outlook
In this work, I have used an orbit-fitting procedure similar to those developed in Gillessen et al. (2017);Boehle et al. (2016) to derive specific constraints on the size of a DM spike for given outer DM halos dynamically constrained by larger scale observations. These limits are the best direct constraints on a DM spike at the GC for non-annihilating DM and exclude a spike with radius greater than a few tens of pc for cuspy outer halos and a few hundred pc for cored outer halos. The addition of the 2017-2018 data from the VLT, which has monitored the pericentre passage of S2 (GRAVITY Collaboration 2018), will make these constraints significantly more stringent, especially thanks to the impressive capabilities of the imaging NACO instrument, the SINFONI spectrometer, and the exquisite astrometric precision of the GRAVITY instrument. However, I postpone the study of the subsequent constraints on the DM profile at the GC to a future work since additional subtleties related to the relativistic effects that have been detected using the new data (GRAVITY Collaboration 2018) may appear, and this warrants a dedicated study. The problem is also complicated further by having to model two pericentre passages when accounting for the entire data set since 1992, which increases the computing time.
Additional improvements on the data could lead to even stronger constraints on the very inner DM profile at the GC. Firstly, in principle, S stars located further out than S2 would be more suited to probe the extended DM distribution for which the mass increases with radius. However, this comes at the price of longer periods, so that unlike S2, no additional stars have been monitored for about 1.5 periods. As a result, our constraints do not improve when including other stars further out such as S1 or S13 for which no significant precession is detectable yet. However, the situation will change when complete orbits are recorded for these stars. In addition, even more accurate astrometric and spectroscopic data will be instrumental to further improve upon these constraints. In particular, a 30 m extremely large telescope (ELT) should be able to probe an extended mass component as low as a few 10 3 M , i.e. about one order of magnitude better than the current sensitivity, after only 10 years of observation (Weinberg et al. 2005). Moreover, an ELT would be able to break the degeneracy between relativistic effects and precession from an extended mass component. This would translate into sensitivity to DM spikes as small as a few pc even for cored outer halos, and even smaller for steeper halos.

Appendix B: Best-fit parameters
Here for completeness I provide the posterior probability distributions obtained for the BH-only model and the BH+spike model with an outer halo of slope γ = 0.25.

Appendix C: Parameters, observables, and analytic solution
Here I recall textbook results of standard mechanics in a central potential (e.g., Bate et al. 1971). The geometry of the problem is described in Fig. C.1. The orbit of the star lies on a plane due to conservation of angular momentum. Let r and θ be the polar coordinates of the star on that plane, and x, y, z the associated Cartesian coordinates. For a stellar orbit we have x = r cos θ, y = r sin θ and z = 0. The x axis is defined in such a way that it points towards the pericentre of the orbit. 7 The dynamics in the orbital plane is well known from first principles. However, observations of the position and velocity of the star are made on the plane of the sky orthogonal to the line of sight pointing to the central object of interest, namely the SMBH Sgr A*. Therefore, we need to transform the phase-space coordinates of the star to a reference coordinate system associated with the plane of the sky. Let X, Y and Z be the Cartesian coordinates in that reference system. Conventionally, the positive X direction points to the North and the positive Y direction to the East, while the positive Z direction points to the Earth. As shown in Fig. C.1, the orbital plane is inclined with respect to the plane of the sky at an angle I. The orbit of the star crosses the plane of the sky in the positive Z direction at the so-called ascending node. The angle between the directions of the ascending node and the pericentre is denoted as ω. Finally, the angle between the X axis-which serves as a reference direction-and the direction of the ascending node is referred to as the longitude of the ascending node, Ω.
As shown in Fig. C.1, the coordinates associated with the plane of the sky are connected to the coordinates in the plane of the orbit via three rotations characterized by ω, I and Ω: first a rotation through −ω around the z axis, then a rotation through I around the axis connecting the central object and the ascending node, and a rotation through −Ω around the Z axis. This is expressed in matrix form as: This can be rewritten as the following functions of time t: where I have used x(t) = r(t) cos θ(t), y(t) = r(t) sin θ(t) and z(t) = 0 for the star. Now, the actual observables are angular coordinates on the sky-more specifically right ascension α and declination δ-as well as radial velocity v r . By definition, right ascension and declination are related to the X and Y coordinates on the sky via α = Y/R 0 and δ = X/R 0 , where R 0 is the distance between Earth and Sgr A*. 8 The motion of the BH is then accounted for through a linear term in the angular position of the star as a function of time. The observable angular coordinates of the star as a function of time (α * (t), δ * (t)) are thus given by where α * /BH (t) = Y(t)/R 0 and δ * /BH (t) = X(t)/R 0 are the angular coordinates of the star in the frame of the BH-with X(t) and Y(t) given by Eqs. (C.2) and (C.3), respectively-and t ref is a reference time, taken to be 2009 yr as in Gillessen et al. (2017). The motion of the BH also causes a shift in radial velocity of the star which reads v r, * ,IR = −Ż * + v r,BH , (C.7) where˙≡ d/dt. 9 10 8 Note that we are in the regime of small angles. 9 Note the minus sign in Eq. (C.7) which accounts for the definition of the Z axis. 10 Measured radial velocities are given in the reference frame of the LSR.