Exact expressions for the pericenter precession caused by some Dark Matter distributions and constraints on them from orbital motions in the Solar System, in the double pulsar and in the Galactic center

We analytically calculate the secular precession of the pericenter of a test particle orbiting a central body surrounded by a continuous distribution of Dark Matter (DM) by using some commonly adopted density profiles for it. We obtain exact expressions without resorting to a-priori simplifying assumptions on the orbital geometry of the test particle. Our formulas allow us to put constraints on the parameters of the DM distributions considered in several local astronomical and astrophysical scenarios, such as the Sun's planetary system, the double pulsar, and the stellar system around the supermassive black hole in Sgr A$^{\ast}$, all characterized by a wide variety of orbital configurations.


Introduction
Latest results [1] from the ESA's Planck mission [2] have yielded a significative revision of the currently accepted content of non-baryonic Dark Matter (DM) of the Universe, which would now amount to about 26.5 ± 1.1% with respect to the previously accepted value of 23.2 ± 1.7% from WMAP data. Indeed, the Planck-based cold DM density, normalized to the critical density, is [1] Ω c h 2 = 0.1199 ± 0.0027, with h = 0.673 ± 0.012, while the analysis of nine years of WMAP data yielded [3] Ω c h 2 = 0.1138 ± 0.0045, with h = 0.70 ± 0.022.
Moving to a galactic scale, the Large Area Telescope (LAT) on the Fermi Gamma Ray Space Telescope spacecraft recently discovered a gamma-ray excess at the Galactic center [4][5][6] which may be due to DM annihilation phenomena. The center of the Milky Way is one of the Galactic regions where a high DM density is expected. As a consequence, if DM was made of self-annihilating constituents into the Standard Model (SM) particles, signals of DM annihilations should come primarily from the Galactic center. The gamma-ray excess signature recorded by the Fermi-LAT instrument [4][5][6] is in the form of a monochromatic gamma-ray line with an energy of about 130 GeV. Later, a second peak at 110 GeV was discovered [7]; the same double peak spectrum was independently observed also in gamma-ray excess from nearby galaxy clusters [8]. The occurrence of a double peak is a generic prediction of DM annihilation pattern in gauge theories.
There is a lingering debate concerning the amount of DM in our Solar System [9][10][11][12]. Recent, controversial estimates by Moni Bidin et al. [10] of the DM local density in the neighbourhood of our Solar System, based on certain assumptions about the vertical velocity dispersion of old tracer stars of the thick disk at the Sun's neighborhood, point towards a value of ρ DM < 0.04 GeV cm −3 . (1) The bound of eq. (1) is smaller by about one order of magnitude than the usually accepted estimate The figure in eq. (2) has been considered as still valid in other recent studies [11,12] criticizing Moni Bidin et al. [10]. Indeed, Bovy and Tremaine [11], who retained the assumptions by Moni Bidin et al. [10] implausible, obtained by using a different hypothesis concerning the mean azimuthal velocity of the stellar tracers, while Nesti and Salucci [12] inferred ρ DM = 0.43 ± 0.11 ± 0.10 GeV cm −3 with a different method, independently of the Galactic DM density profile. An accurate knowledge of the DM distribution in the neighbourhood of the Solar System is relevant for the attempts aimed to directly detect DM particles in laboratory-based experiments such as, e.g., CDMSI [13], CDMSII [15], DAMA/NaI [16] and its successor DAMA/LIBRA [17], XENON10 [18] and ZEPLIN III [19].
Interestingly, the silicon (Si) detectors of the CDMS II experiment has recently revealed three WIMPcandidate events at a 3σ level. [14]. Moving to the Earth's neighbourhood, hints of DM might have been detected [20] by the Alpha Magnetic Spectrometer (AMS) on the International Space Station, which has recently measured an anomalous high-energy positron excess in Earth-bound cosmic rays [21]. Thus, it has become even more important to devise independent means to effectively gain information on the constitution and the distribution of such a hypothesized, elusive ingredient of the natural world. In Section 2, we will analytically calculate the secular perturbations to the motion of a test particle orbiting a central body surrounded by a continuous DM distribution. We will adopt some spherically symmetric local DM density profiles [22] by computing the corresponding corrections U DM to the Newtonian twobody potential through the Poisson equation. Then, we will evaluate them onto the unperturbed Keplerian ellipse of the test particle, and we will average them over one orbital period P b of the test particle in order to use the standard Lagrange planetary equations [23] for the variation of the Keplerian orbital elements. We will, thus, be able to straightforwardly infer exact expressions of the orbital rates of change due to the different DM distributions considered without resorting to any a-priori simplifying assumptions on the orbital configuration of the test particle. Such formulas will be useful to constrain the DM local distributions in different astrophysical and astronomical scenarios such as planetary systems, binaries hosting compact objects, stellar systems around supermassive black holes, etc. characterized by a wide variety of orbital geometries. Only more or less approximate expressions for the DM-induced pericenter precession exist in literature.
In Section 3 we will apply them to the planets of the Solar System, to the double pulsar, and to the stars at the Galactic center in Sgr A * . Section 4 summarizes our findings.
In many cases, more or less approximate expressions for the anomalous perihelion precession induced by certain spherically symmetric DM distributions were used, in particular by considering nearly circular orbits. In this Section, we will overcome such a restriction by calculating exact expressions which can, thus, yield more accurate constraints in view of the increasing level of accuracy in determining the orbits of some of the major bodies of the Solar System. Moreover, our results can be used also with systems characterized by highly eccentric orbits such as extrasolar planets, binaries hosting compact objects, and stellar systems orbiting supermassive black holes lurking in galactic nuclei.

Exponential DM density profile
By choosing an exponentially decreasing DM profile [36] ρ DM (r) = ρ 0 exp − r λ , where ρ 0 is the density at a distance r equal to the characteristic scale length λ, the Poisson equation yields The density profile of eq. (5) is a particular case of the Einasto profile [38], often adopted to describe DM halos in galaxies [22]. The effect of eq. (6) on the orbital motion of a test particle moving around a central body of mass M surrounded by a DM continuous distribution characterized by the density profile of eq. (5) can be computed perturbatively by assuming eq. (6) as a small correction to the Newtonian monopole U N (r) = −GM/r. The average of eq. (6) over one orbital revolution of the test particle turns out to be In eq. (7), I 0 (x) , I 1 (x) are the modified Bessel functions of the first kind, and x . = ae/λ, where a, e are the semimajor axis and the eccentricity of the orbit of the test particle, respectively. The Lagrange equation for the variation of the longitude of the pericenter ̟ [23] and eq. (7) yield the following secular precession where n b . = GM/a 3 = 2π/P b is the Keplerian mean motion of the test particle. Although Pitjev and Pitjeva [36] considered the density profile of eq. (5), they did not explicitly calculate the resulting pericenter precession. The precession of eq. (8) is an exact result in the sense that neither a-priori simplifying assumptions on e nor on λ were assumed. Note that eq. (8) is not singular in e since, in the limit e → 0, it reduces to In the case of an infinite length scale, i.e. for an uniform mass density, eq. (8) reduces to the known result [27][28][29]36]

Power-law DM density profile
By assuming a power-law DM distribution profile usually used for the galactic halos and in several DM-related studies [22,39], the Poisson equation yields The orbital effects of eq. (11) have been computed more or less explicitly and at various levels of approximations in [30,32,34,35,40,41]. Actually, also the averaged pericenter precession induced by eq. (12) can be exactly computed with the Lagrange pertubative scheme [23] without resorting to any simplifying assumptions on the eccentricity e of the test particle.
It turns out to be where the cumbersome expressions of the coefficients p j (e; γ, λ) , j = 1, 2, . . . 16 are explicitly displayed in Appendix A. It is worthwhile noticing that eq. (13) is not any sort of somewhat truncated series expansion, and that the coefficients in Appendix A are exact in the sense that no a priori simplifying assumptions were assumed in calculating them. By performing a series expansion of eq. (13) in powers of e one obtains The first term in eq. (14) yields a perihelion shift per orbit in agreement with that by Frère et al. [32].

Confrontation with the observations
In this Section, we will use our analytical predictions of eq. (8) and eq. (13) to infer constraints on ρ 0 as functions of λ and γ in some astronomical and astrophysical scenarios.

Planets of the Solar System
By following the approach adopted in several researches on DM [26][27][28][29][30][31][32][33][34][35][36] and on other putative nonstandard dynamical effects [27,42] in our Solar System, we will confront our theoretically predicted pericenter precessions of eq. (8) and eq. (13) with the latest determinations of the admissible ranges ∆̟ for any possible anomalous perihelion precessions obtained by fitting up-to-date dynamical models to a centennial record of planetary observations of several types. The term "anomalous" refers to the standard Newtonian-Einsteinian dynamics, fully modeled in the most recent planetary dynamical theories. At present, two independent teams of astronomers are engaged in producing, among other things, such corrections ∆̟ to the planetary perihelion precessions [36,43]; here we will adopt the latest estimates by Pitjev and Pitjeva [36], summarized in Table 1, which are based on the most recent version of the EMP ephemerides. DM was not explicitly modeled in the EPM routines; thus, the supplementary precessions of Table 1 are well suited to put constraints on it in a phenomenological, model-independent way. Here and in the rest of the paper, "model-independent" has to be intended as independent of the nature and of the physical properties of the DM particles. On the other hand, the constraints inferred in this paper Table 1. Supplementary precessions ∆̟ of the longitudes of the perihelion for some planets of the Solar System estimated by Pitjev and Pitjeva [36] with the EPM2011 ephemerides. They processed 676804 data of various kinds covering about one century (1913-2010) by estimating more than 260 parameters. While some tracking data from Cassini were used, those from Messenger, currently orbiting Mercury, were not included in the data record. Pitjev and Pitjeva [36] fully modeled all standard Newtonian-Einsteinian dynamics, apart from the general relativistic Lense-Thirring effect caused by the Sun's rotation. However, its expected magnitude from general relativity is smaller than the uncertainties quoted here for the planets used in the text (Earth, Mars, Saturn), so that it does not impact our results. . The units are milliarcseconds per century (mas cty −1 ). The errors released σ ∆̟ , in mas cty −1 as well, generally exceed the formal, statistical ones by several times [36]. depend on how DM clusters in space, i.e. on the specific density profiles adopted. Moreover, by assuming that the whole ranges ∆̟ in Table 1 are entirely due to DM, the resulting constraints will be relatively "generous", mitigating the risk of getting unrealistically tight bounds. A complementary approach which could be followed, recommended by some researchers [44] when non-standard effects such as DM are involved, consists in explicitly modeling DM and fitting such ad-hoc modified dynamical theories to the same data set in order to estimate a dedicated solved-for parameter in a least-square way. One could try with different DM density profiles to establish the one providing the best fit to the data. Such a strategy would be particularly meaningful in the case where non-zero extra-precessions, at a statistically significant level, should be extracted from future observations.

Constraints for the exponential density profile
Here, we deal with the exponential density profile of eq. (5).

Constraints for the power-law density profile
Here, the power-law distribution of eq. (11) is considered.  Table 1 for Earth, Mars and Saturn were adopted along with the analytical prediction of eq. (8). The characteristic length λ, in AU, is assumed to vary from r min = a(1 − e) to r max = a(1 + e).  Figure 2 displays the bounds on ρ 0 , in GeV cm −3 , as a function of the characteristic scale length λ and of the parameter γ for the Earth, Mars and Saturn. Also in this case, λ covers the orbit extensions of the planets considered, while [32,35] 0 ≤ γ ≤ 4. Also in this case, Saturn yields the tightest bounds: 5 × 10 3 GeV cm −3 (γ = 0) ≤ ρ 0 ≤ 8 × 10 3 GeV cm −3 (γ = 4) corresponding to 8.9 × 10 −21 g cm −3 ≤ ρ 0 ≤ 1.4 × 10 −20 g cm −3 . They are about one order of magnitude better than those inferred in Section 3.1.1. .

The double pulsar
Thanks to their high matter density, astrophysical compact objects such as neutron stars should be able to efficiently capture significative amounts of DM [45][46][47][48][49][50][51][52][53][54][55] depending on the nature of its particles, on how it clusters in space and on the history of the stars themselves. Concentrations of DM around neutron stars may give rise to phenomena such as gamma-ray emission [54,55] which, on the one hand, depend on the physical properties of DM particles themselves [56] and, on the other hand, can be generated also by other physical mechanisms [57,58]. Thus, it is valuable to constrain the DM density around binaries hosting neutron stars in a phenomenological and model-independent way.

The stellar system orbiting the Galactic black hole in Sgr A *
The issue of the DM distribution at the center of galaxies in presence of supermassive black holes has been treated in, e.g., [40,41,[63][64][65][66]. It depends on several factors such as [63] the galactic halo density profile itself [67] and the nature of the DM particles themselves. Moreover, as pointed out in [66], significant differences in the final DM distribution close to the black hole are found depending on the theoretical scheme adopted for the calculation.
So far, the presence of DM at the center of the Milky Way has been indirectly guessed from certain phenomena such as gamma-ray emissions interpreted as DM annihilation [68,69]. However, such interpretations are not free from more or less tested assumptions about several other concurring physical phenomena [70,71]. Recently, the Fermi satellite detected a gamma-ray excess at the Galactic center which may be due to DM annihilation phenomena [4][5][6]. Thus, it is important to constrain DM distributions at the Galactic center in a dynamical, model-independent way.  Table 1 for Earth, Mars and Saturn were adopted along with the analytical prediction of eq. (13). The characteristic length λ, in AU, is assumed to vary from r min = a(1 − e) to r max = a(1 + e).   To this aim, we will look at the orbital motions of the stellar system [39,72,73] revolving about the supermassive black hole located at the Galactic center in Sgr A * [74]. So far, available data cover one full orbital revolution for two stars: S0-2, having an orbital period P b = 16.17 ± 0.15 yr and eccentricity e = 0.898 ± 0.005 [73], and S0-102, characterized by P b = 11.5 ± 0.3 yr and eccentricity e = 0.68 ± 0.02 [73]. The mass of the black hole is M • = (4.1 ± 0.4) × 10 6 M ⊙ [73]. From such figures, a naive, order-of-magnitude evaluation on the accuracy which could be reached in determining the stellar perinigricon precessions can be made; for S0-2 we have σω ≈ 0.6 deg yr −1 , while for S0-102 it can be inferred σω ≈ 4.8 deg yr −1 . Then, we will adopt S0-2 as probe to constrain DM density as done in the previous Sections. In view of its high eccentricity, the exact formulas of eq. (8) and eq. (13) are well suited for such a task. From Figure 4 it turns out that the bounds for the exponential density vary sensibly over the S0-2 orbit extension in view of its high eccentricity; indeed, it is ρ 0 ≈ 4 × 10 15 GeV cm −3 , corresponding to 7 × 10 −9 g cm −3 , for λ = r min , while it reduces to ρ 0 ≈ 2 × 10 13 GeV cm −3 , corresponding to 3 × 10 −11 g cm −3 , for λ = r max . As pointed out in [39], as far as supposedly baryonic DM is concerned, both theoretical [75,76] and observational [77] motivations for considering a power-law profile such as eq. (11) at the galactic centers exist in literature. The same kind of potential was considered for non-baryonic DM as well [41,64]. According to Figure 4, it yields 1.2 × 10 13 GeV cm −3 (γ = 0) ≤ ρ 0 ≤ 1 × 10 16 (γ = 4, λ = r min ) GeV cm −3 corresponding to 2.1 × 10 −11 g cm −3 ≤ ρ 0 ≤ 1.8 × 10 −8 g cm −3 .
For an earlier idea of constraining the DM concentration near the Galactic center from the perinigricon of S0-2, see [40,41].

Summary and conclusions
Recent results concerning the amount of non-baryonic DM, both at global and local scales, stress the need of further deepening the research of accurate and independent strategies to gain information about the distribution of such a hypothetical key ingredient of the natural world.
To this aim, we looked at the effects induced by some spherically symmetric DM density profiles on the motion of a test particle orbiting a localized body surrounded by an extended DM distribution. In view of the increasing accuracy in the determination of the orbits of the planets of our Solar System and of the possibility of looking also at different astronomical and astrophysical laboratories characterized by a wide variety of orbital configurations, we analytically calculated the DM-induced pericenter precession without resorting to any a-priori simplifying assumptions concerning the orbital geometry of the test particle. In this respect, our results, obtained perturbatively by means of the Lagrange planetary equations, are exact, being valid for any value of the eccentricity of the orbit of the test particle.
Our results can be used in future when new, more accurate data will be collected and processed. As a complementary approach which could be followed, DM dynamical effects should be explicitly modeled in the softwares used to reduce the planetary observations, and dedicated solved-for parameters should be estimated in fitting the newly constructed models of the forces acting on the planets to given data records.