BBinary pulsars as dark-matter probes
Paolo Pani
1, 2, ∗ Dipartimento di Fisica, “Sapienza” Universit`a di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 5, 00185, Roma, Italy. CENTRA, Departamento de F´ısica, Instituto Superior T´ecnico,Universidade de Lisboa, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal.
During the motion of a binary pulsar around the Galactic center, the pulsar and its companionexperience a wind of dark-matter particles that can affect the orbital motion through dynamicalfriction. We show that this effect produces a characteristic seasonal modulation of the orbit andcauses a secular change of the orbital period whose magnitude can be well within the astonishingprecision of various binary-pulsar observations. Our analysis is valid for binary systems with orbitalperiod longer than a day. By comparing this effect with pulsar-timing measurements, it is possibleto derive model-independent upper bounds on the dark-matter density at different distances D fromthe Galactic center. For example, the precision timing of J1713+0747 imposes ρ DM (cid:46) GeV / cm at D ≈ D (cid:46)
10 pc could provide stringent constraints ondark-matter halo profiles and on growth models of the central black hole. The Square KilometerArray can improve current bounds by 2 orders of magnitude, potentially constraining the localdensity of dark matter to unprecedented levels.
PACS numbers: 95.35.+d, 97.60.Gb, 95.30.Cq.
I. INTRODUCTION
Dark matter (DM) is five times as abundant as bary-onic matter in the Universe but its properties remainmysterious. While there is strong evidence for DM par-ticles to be nonrelativistic, their mass, spin, charges,and annihilation and interaction cross sections are un-known [1]. The plethora of DM candidates makes directdetection extremely challenging and model dependent. Afurther astrophysical problem related to DM is to mea-sure its density profile in the Milky Way. Conventionalcold DM cosmological models predict cuspy halos whichmight give rise to strong gamma-ray emission throughDM annihilation. However, such DM cusps seem to bein conflict with observations (see e.g. Refs. [2, 3]).In this intricate scenario, the universality of the grav-itational interactions might represent a stronghold formodel-independent tests of DM. In this paper we pointout that binary pulsars are the ideal laboratory for suchgravitational tests. Pulsar-timing techniques allow us tomeasure some of the orbital parameters with exquisiteprecision, and are routinely used to measure the mass ofneutron stars (e.g. Refs. [4, 5]) and to perform some ofthe most stringent tests on Einstein’s theory of generalrelativity, as in the case of the celebrated Hulse-Taylorbinary pulsar [6] and of the more recent double-pulsarsystem [7] (see also Ref. [8] and Refs. [9–11] for somereviews).During the motion of a binary pulsar around the Galac-tic center, the pulsar and its companion experience awind of DM particles similar to the one that producesan annual modulation in the scattering rate in direct-detection experiments due to Earth’s motion around the ∗ [email protected] Sun [12, 13]. We argue that the effects of DM dynamicalfriction (i.e., the drag force due to the gravitational in-teraction of the orbiting bodies with their wakes [14–16])on the binary motion can produce observable effects andcan therefore be used to put stringent constraints on theDM density in our galaxy.
II. DM EFFECTS ON THE TWO-BODYPROBLEM
We consider a binary system with masses m and m whose center of mass moves through a DM distributionwith a roughly constant local density, ρ DM . In the center-of-mass frame, and neglecting relativistic effects for sim-plicity, the equations of motion for the two-body systemread m i ¨ r i = ± Gm m r r + F DF i , (1)where the upper (lower) sign refers to i = 1 ( i = 2), r i is the position vector of the mass m i , r := r − r is therelative position of the bodies, and F DF i is the gravita-tional drag on the i -th object. The effects of DM accre-tion are some orders of magnitude smaller than the effectdiscussed here and will be neglected.Dynamical friction depends on the nature of themedium and on its gravitational interaction with the ob-jects. In the case of DM, the medium can be modeledas a collisionless gas as long as the DM mean free pathis much larger than the size R i of the objects. This as-sumption requires σ DM m DM (cid:28) (cid:18) R (cid:12) R i (cid:19) (cid:18) GeV / cm ρ DM (cid:19) cm / g , (2)where σ DM and m DM are the DM self-interaction crosssection and the mass of the DM particles. DM self-interactions are constrained by several observations, a r X i v : . [ a s t r o - ph . H E ] D ec FIG. 1. An illustration of the binary system under consider-ation in the Galactic reference frame. The center of mass ofthe binary moves through a constant-density DM distribution(hatched area) with velocity v w , experiencing a wind of DMparticles with velocity − v w . a conservative bound being σ DM /m DM (cid:46) cm / g (cf.Ref. [17] and references therein). Therefore, even withthe extreme values adopted in Eq. (2), DM is perfectlycollisionless.For a single object in linear motion through a homo-geneous collisionless medium, dynamical friction is gov-erned by Chandrasekhar’s formula [14–16] F DF i = − πρ DM λ G m i ˜ v i (cid:18) erf ( x i ) − x i √ π e − x i (cid:19) ˜ v i , (3)where ˜ v i = ˙ r i + v w is the velocity of the i -th body relativeto the DM wind, x i := ˜ v i / ( √ σ ), σ is the dispersion ofthe DM Maxwellian velocity distribution, and λ ≈ +10 − is the Coulomb logarithm [15, 16]. We have also intro-duced the rotational velocity of the binary around thegalaxy, v w = v w (cos α sin β, sin α sin β, cos β ), which (ne-glecting the rotational velocity of the DM halo) is oppo-site to the velocity of the wind of DM particles relativeto the center of mass (cf. Fig. 1 for an illustration of thesystem under consideration). Over the time scale of ob-servations (at most some decades) v w is approximatelyconstant.Although Eq. (3) was derived for linear motion, ityields remarkably precise results also for the case of ageneric orbit [18–20]. When applied to a binary sys-tem, Chandrasekhar’s formula neglects the interactionof one component with its companion’s wake. This ef-fect is negligible when the characteristic size L of thewake is smaller than the orbital radius r and when theorbital velocity v is sufficiently small as to give time tothe wake to disperse after an orbit. By estimating L from the size of the gravitational sphere of influence, mσ ∼ Gm i m DM /L , the first condition L (cid:28) r reads [21] P b (cid:29) Gm i σ ∼ . (cid:18) m i . M (cid:12) (cid:19) (cid:18)
150 km / s σ (cid:19) day , (4) where P b is the orbital period of the binary. This con-dition also implies v (cid:28) σ , which guarantees that theparticles in the wake of one object disperse before thearrival of the companion along the orbit. Thus, at leastfor binaries at large orbital distance, the superposition ofEq. (3) for i = 1 , P b ∼ . σ (cid:38)
200 km / s. With these motivations, here we re-strict to the condition (4) and adopt Eq. (3) to modelthe gravitational drag. A more general analysis valid forany orbital period is underway. Some limitations andthe range of validity of Eq. (3) in binary systems arediscussed in Appendix A. A. Perturbed orbits
By introducing the center-of-mass position vector R :=( m r + m r ) /M (with M = m + m ), the system (1)( i = 1 ,
2) can be conveniently written as a set of twoequations of motions for r and R . By rewriting Eq. (3)as F DF i = − Ab i m i M ˜ v i , where A := 4 πρ DM λG M and b i := 1˜ v i (cid:20) erf ( x i ) − x i √ π e − x i (cid:21) , (5)it is straightforward to show that the equations of motionreduce to ˙ v = − GMr r + a η v + a ( v w + V ) , (6)˙ V = a η v + a ( v w + V ) , (7)where we have defined v = ˙ r , V = ˙ R , a = − A ( b + b ) , (8) a = A b ∆ + + b ∆ − ] , (9) a = − A (cid:2) b ∆ + b ∆ − (cid:3) , (10)∆ ± = ∆ ±
1, ∆ = √ − η , η = µ/M , and µ = m m /M is the reduced mass. Note that, because of the externalforce F DF i , the center of mass is accelerated, ¨ R (cid:54) = 0.The dynamics can be greatly simplified by treatingthe contributions of dynamical friction perturbatively.This assumption is perfectly justified by the tiny magni-tude of these corrections relative to the Keplerian terms.To compute the first-order (in the DM density) correc-tions, we adopt the formalism of the osculating orbits (cf.Ref. [22] for a modern treatment). For simplicity, we fo-cus on the corrections to a circular binary; the ellipticalcase is a straightforward generalization.To zeroth order, the motion is planar and it isconvenient to adopt polar coordinates where r = r ( t )(1 , ϕ ( t ) ,
0) and Ω = ˙ ϕ is the orbital angular veloc-ity. Obviously, R is constant and the motion is Keplerian, r Ω = ( GM Ω ) / = v .To compute the first-order corrections, all quantitiesthat are multiplied by ρ DM can be evaluated using theKeplerian solution. In particular, to zeroth order V = V = 0 and the right-hand side of Eqs. (6) and (7) canbe further simplified. To first order, Eq. (6) has the form˙ v = − GMr r + f , with f = a η v + a v w , and does notdepend on the motion of the center of mass. FollowingRef. [22], the adiabatic-evolution equations for the orbitalparameters in the quasicircular case are˙ a = 2 (cid:114) r GM S ( t ) , (11)˙ e = (cid:114) r GM [ R ( t ) sin Ω t + 2 S ( t ) cos Ω t ] , (12)˙ ι = (cid:114) r GM W ( t ) cos(Ω t + ω ) , (13)˙Ω = 1sin ι (cid:114) r GM W ( t ) sin(Ω t + ω ) , (14)where a , e , ω , ι and Ω are the semiaxis major, the eccen-tricity, the longitude of pericenter, the inclination, andthe longitude of the ascending node, respectively. Thesource terms in the equations above read R := f · n = a n · v w = a v w sin β cos(Ω t − α ) , (15) S := f · λ = a η λ · v + a λ · v w = a ηv − a v w sin β sin(Ω t − α ) , (16) W := f · e z = a v w cos β , (17)where we used the definitions of the orbital basis vectorsto zeroth order, namely n = (cos Ω t, sin Ω t, λ =( − sin Ω t, cos Ω t, e z = (0 , , P b := 2 π/ Ω .A relevant quantity is the time derivative of the orbitalperiod which – in a quasiadiabatic approximation valid inthe perturbative regime – can be obtained from Eq. (11)and from Kepler’s third law, ˙ P b = P b r ˙ a . Using Eq. (11),we obtain˙ P DF b ( t ) = 3 vP b a η − a Γ sin β sin(Ω t − α )] , (18)where Γ := v w /v . We focus on the secular changes ofthe orbital parameters, so henceforth we shall considerorbital-averaged quantities obtained from Eqs. (11)–(14),e.g. (cid:104) X (cid:105) := P − b (cid:82) P b dtX ( t ). B. Other contributions to ˙ P b In addition to the changes of the orbital period in-duced by the secular evolution of the semiaxis major [cf.Eq. (18)], the observed variation of P b in a binary systemis affected by kinematic effects, namely by the center-of-mass acceleration along the line of sight and by thevariation of the orbital inclination ι . Both effects pro-duce a Doppler shift of pulse frequencies which affectsthe measurement of P b (cf., e.g., Ref. [23]). The induced change in the orbital period due to thecenter-of-mass acceleration reads˙ P cm b = P b ˙ V · e Z , (19)where e Z is the unit vector parallel to the line of sightas defined in Ref. [22]. In Appendix B we compute thecontribution above, whose final expression reads˙ P cm b ( t )= vP b { a η sin ι cos(Ω t + ω )+ a Γ [cos ι cos β + sin ι sin β sin( α + ω )] } . (20)Finally, the contribution of the variation of the or-bital inclination reads ˙ P ιb = tan ιP b ˙ ι which, by usingEqs. (13) and (17), reduces to˙ P ιb ( t ) = 32 a Γ P b tan ι cos β cos(Ω t + ω ) . (21)Note that ˙ P ιb ( t ) = 0 when ι = 0 (i.e. when the orbit isperpendicular to the line of sight) and when β = π/ P DF b +˙ P cm b + ˙ P ιb , obtained from Eqs. (18), (20) and (21). In ourcase the contribution of ˙ P cm b is typically subdominantand can be neglected, whereas ˙ P ιb can be of the sameorder of ˙ P DF b . C. Analytical results
In generic situations the orbital average of Eqs. (18),(20), and (21) should be performed numerically becausethe quantity b i ( t ) defined in Eq. (5) (and evaluated usingthe Keplerian solution) is a cumbersome function of thetime. Before presenting the results of such numericalintegration in Sec. II D, here we discuss some limits inwhich the secular changes can be computed analytically.For simplicity we focus on the ˙ P DF b term; the corrections˙ P cm b and ˙ P ιb can be computed using the same procedure.
1. Large- σ limit In the limit σ (cid:29) v w , v , from Eq. (5) we obtain b i ∼ (cid:114) π σ × (cid:26) − v σ (cid:2) + ∆ ∓ − ∓ sin β sin(Ω t − α ) (cid:3) + . . . (cid:27) , where we have used the Keplerian solution to lowest or-der, and the upper (lower) sign refers to i = 1 ( i = 2).Note that, to leading order, b ∼ b ∼ const. By insertingthe above expression into Eq. (11), it is straightforwardto obtain the analytical result˙ P DF b ( t ) ∼ − √ πG µλρ DM P b σ (cid:26) η sin β sin(Ω t − α )++ 3 v ησ (cid:2) Γ η cos(2Ω t − α ) + Γ∆ (cid:0) η − Γ (cid:1) sin β sin(Ω t − α ) + η (cid:0) η − − + 2Γ cos 2 β sin (Ω t − α ) (cid:1)(cid:3)(cid:27) , where the terms on the last line are the next-to-leading-order corrections in the large- σ limit. The orbital changeas a function of time has a characteristic modulation withperiod P b , an amplitude proportional to the DM windvelocity, and a phase given by the angle α . This modula-tion is expected, due to the periodic motion at Keplerianorder. However, to leading order the oscillatory term av-erages to zero over an orbital period. In this limit thesecular decay, (cid:104) ˙ P DF b (cid:105)∼ − √ πG µλρ DM P b σ × (cid:20) v σ (2 η − − + Γ cos 2 β ) (cid:21) , (22)does not depend on the DM wind v w to leading order.Equation (22) extends the result by Gould [24] obtainedin the large- σ regime within a rather different framework;in fact, the leading-order term of the expression above matches exactly that derived in Ref. [24]. This agreementgives further confirmation of the validity of Eq. (3), atleast in the σ (cid:29) v limit.Note also that, to leading order, a in Eq. (21) is aconstant and therefore (cid:104) ˙ P ιb (cid:105) = 0. Likewise, when β = 0the orbital average (cid:104) ˙ P ιb (cid:105) = 0 in general (and not onlyin the large- σ limit), because in such case a does notdepend on time.
2. Large- v w limit Another relevant limit is v w (cid:29) σ, v , which yields b i ∼ v w (cid:20) ∓
2Γ sin β sin(Ω t − α ) (cid:21) , (23)where again the upper (lower) sign refers to i = 1 ( i = 2).From this expression it is straightforward to obtain˙ P DF b ( t ) ∼ − G M π ∆ λρ DM vv w Ω (cid:34) sin β sin(Ω t − α ) + 2 η (cid:2) − β sin (Ω t − α ) (cid:3) Γ∆ + 3 η sin β sin(Ω t − α )Γ (cid:35) . The O (1 /v w ) and O (1 /v w ) terms are zero when β = 0and their orbital average vanishes for any β . Therefore,in this limit the secular change of the orbital period reads (cid:104) ˙ P DF b (cid:105) ∼ − πG M P b ηλρ DM (1 + 3 cos 2 β ) v w . (24)The leading-order coefficient depends on the direction ofthe DM wind but not on the DM velocity dispersion.Interestingly, (cid:104) ˙ P DF b (cid:105) becomes positive when 55 ◦ (cid:46) β (cid:46) ◦ . In this case the orbit tends to outspiral (i.e., toget softer [25]) due to the effect of DM dynamical fric-tion. This configuration seems to violate Heggie’s law [26]according to which hard binaries (i.e. those with bind-ing energy E b (cid:29) m DM σ ) get harder ( (cid:104) ˙ P DF b (cid:105) <
0) in amedium of fast low-mass objects. This result might haveimplications in studies of Galactic dynamics and we planto investigate it more in detail in a separate publication.Our result is not in contrast with the findings of Ref. [24],because the latter were obtained for v w = 0, whereas herewe assume v w (cid:29) σ, v . Note that, in addition to the large- v w limit, the condition (4) also needs to be satisfied.
3. Small- v w limit In the limit v w ∼
0, we obtain˙ P DF b ( t ) ∼ − GM µλρ DM P b √ πm m σ × (cid:8) M √ π (cid:2) m erf( y ) + m erf( y ) (cid:3) σ −√ m m (cid:16) e − y m + e − y m (cid:17) v (cid:111) . (25)where y i := m i v √ Mσ . Since in this case ˙ P DF b does notdepend on time, the orbital average coincides with theequation above, which also reduces to the leading-orderterm of Eq. (22) when σ (cid:29) v . Note that in this case (cid:104) ˙ P DF b (cid:105) <
0, in agreement with the results of Ref. [24] andwith Heggie’s law.
D. Numerical results
In various situations of phenomenological interest v ∼ v w ∼ σ ∼ O (100 km / s) and the analytical limits dis-cussed above are not valid. In this case one needs toresort to a numerical integration, which is in any case -16 -15 -14 -13 -12 | P b . D F | σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/sanalytical w [km/s]10 -16 -15 -14 -13 - P b . D F β = π /2 β =0 FIG. 2. Secular change of the orbital period of a binary pulsardue to DM dynamical friction as a function of the DM windvelocity v w , for different values of the DM velocity dispersion,and for β = π/ β = 0 (bottom panel).We assumed the following reference values: P b = 100 day, m = 1 . M (cid:12) , m = 0 . M (cid:12) and ρ DM = 2 × GeV / cm .The analytical result refers to Eq. (24), whereas the horizontallines refer to Eq. (25). straightforward given the simple form of Eqs. (11)–(14).Henceforth the notation (cid:104) X (cid:105) for the orbital average of aquantity X ( t ) will be omitted.In Fig. 2 we present some exact numerical results for˙ P DF b and compare them with Eqs. (24) and (25). Asexpected, Eq. (24) is an excellent approximation to theexact result when v w (cid:29) σ, v . Furthermore, we find thatthe analytical result (25) is accurate in the entire range v w (cid:46) σ and not only when v w →
0, at least for longorbital periods. Note also that ˙ P b can become positiveabove a critical value of v w when β is sufficiently large.This is in agreement with Eqs. (25) and (24). While theformer equation predicts ˙ P b < v w ∼
0, the latterpredicts ˙ P b > ◦ (cid:46) β (cid:46) ◦ and at large valuesof v w . By continuity, there exists a value of v w at which˙ P b changes sign.The secular change ˙ P DF b as a function of the orbital pe-riod P b is shown in Fig. 3 for various choices of the angle β , of the DM wind velocity v w , and of the DM velocitydispersion σ . As expected from Eq. (3), corrections arelarger for longer orbital periods, although in some partic-ular cases (when β ∼ π/
2) the behavior of the function˙ P b might be nonmonotonic An extensive numerical search shows that the secu- As previously discussed, ˙ P b can change sign for large values of β .In the logarithmic scale adopted for the vertical axis of Fig. 3(and of subsequent figures), the zero crossing corresponds to“cusps” which occasionally appear in the function ˙ P DF b when β ≈ π/ lar change ˙ P DF b is independent of α (in agreement withthe analytical results previously presented) and that thecorrections depend only mildly on the value of β . InFigs. 2 and 3 we consider a typical neutron star-whitedwarf system with m ≈ . M (cid:12) and m ≈ . M (cid:12) . Thisis a conservative assumption, since the secular changesgrow with the masses of the binary, as shown in Fig. 4(to be compared with the corresponding Fig. 3).Finally, our perturbative solution is valid for t (cid:28) P b / ˙ P DF b , i.e. for t (cid:28) P b or longer (cf. Figs. 3 and 4)and it is therefore perfectly reliable during the observa-tion time (lasting at most some decades). III. CONSTRAINTS ON THE DM DENSITY
The orbital period of various binary pulsars is found tobe constant over several years within astonishing preci-sion. Such systems are therefore ideal candidates to con-strain the small effect of DM on the binary motion. Forsome systems, measurements of the small secular changesof the orbital period ( ˙ P b , ¨ P b , ...) are also available.The apparent orbital change ˙ P obs b needs to be cor-rected to take into account the differential accelerationin the Galactic gravitational potential [27] and the trans-verse motion of the binary relative to Earth (Shklovskii’seffect [28]). After subtracting these kinematic contri-butions, the intrinsic orbital decay, ˙ P int b , can still becaused by gravitational-wave emission, mass loss, tidaltorques and gravitational quadrupole coupling [29–31](cf. Ref. [23] for an overview).For some binary systems these corrections can be com-puted precisely. For relativistic binaries, the dominantcontribution is typically due to the gravitational-waveemission, ˙ P GW b , and the “excess” orbital decay, ˙ P xs b =˙ P int b − ˙ P GW b , is found to be consistent with zero withinerrors. In other words, for some systems the emission ofgravitational waves as predicted by the quadrupole for-mula of general relativity [32] completely accounts for theobserved intrinsic secular change of ˙ P b . This is the caseof J1738+0333, whose measured excess orbital period, | ˙ P xs b | (cid:46) × − , provides one of the most stringentconstraints on modified theories of gravity [8].Interestingly, the typical amplitude of ˙ P DF b shown inFigs. 2–4 is comparable to or even larger than the ob-served value of ˙ P xs b for some systems [8, 29–31, 33]. As anorder-of-magnitude estimate, the leading term of Eq. (22)can be rewritten as˙ P DF b ≈ − × − µ λ ρ GCDM P (100) b σ (26)where µ = µ/M (cid:12) , λ = λ/ P (100) b = P b / (100 day), σ = σ/ (150 km / s), and ρ GCDM = ρ DM / (2 × GeV / cm )is normalized by approximately the value of the DM den-sity at a distance D ≈ D ≈ -16 -15 -14 -13 -12 - P b . D F σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 -17 -16 -15 -14 -13 - P b . D F β = π /2 β =0v w =100 km/sv w =100 km/s -16 -15 -14 -13 -12 | P b . D F | σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 -17 -16 -15 -14 -13 - P b . D F β = π /2 β =0v w =220 km/sv w =220 km/s FIG. 3. Secular change of the orbital period of a binary pulsar due to DM dynamical friction as a function of the orbitalperiod P b , for different values of the DM velocity dispersion and of the DM wind (left panels: v w = 100 km / s; right panels: v w = 220 km / s), and for β = π/ β = 0 (bottom panels). We considered the same reference values as inFig. 2. -15 -14 -13 -12 | P b . D F | σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 -16 -15 -14 -13 -12 - P b . D F β = π /2 β =0v w =100 km/sv w =100 km/s -16 -15 -14 -13 -12 | P b . D F | σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 -16 -15 -14 -13 - P b . D F β = π /2 β =0v w =220 km/sv w =220 km/s FIG. 4. Same as in Fig. 3 but for m = 2 M (cid:12) and m = 5 M (cid:12) . larger by a factor 100 in the case of adiabatic growth ofthe central BH [35, 36]. Therefore, a binary pulsar with P b ∼
100 day near the Galactic center will display a sec-ular change of the orbital period well within the currentaccuracy of pulsar timing.On the other hand, Eq. (26) can be inverted yielding anapproximate upper bound on the DM density allowed bythe timing measurement of a binary pulsar with orbitalperiod P b and excess orbital change ˙ P xs b , namely ρ DM (cid:46) (cid:32) ˙ P xs b − (cid:33) σ µ λ P (100) b GeV / cm . (27)The bound above is purely indicative, since it is ob-tained using the analytical result (22) valid only in thelimit σ (cid:29) v w , v . Figure 5 shows more precise upper bounds on ρ DM derived from the numerical results by im-posing | ˙ P DF b | (cid:46) × − , i.e. assuming a measurement of˙ P xs b as precise as that for J1738+0333 [8]. Note that thevalues σ = 150 km / s, v w = 220 km / s and σ = 50 km / s, v w = 100 km / s considered in Fig. 5 are representativeof the DM velocity dispersions and Galactic rotation ve-locities near the Solar System and close to the Galacticcenter, respectively.Due to a selection bias, most binary pulsars are de-tected near Earth. Indeed, we searched in the ATNFPulsar Catalogue [37, 38] for those binaries whose (appar-ent) orbital change ˙ P obs b has been measured and found 32systems which are located at a distance D ∈ [1 −
12] kpcfrom the Galactic center. Therefore, even though forsome of these binaries the intrinsic orbital change can be ρ D Mm a x [ G e V / c m ] σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 ρ D Mm a x [ G e V / c m ] β = π /4 β = π /4v w =220 km/sv w =100 km/s ρ D Mm a x [ G e V / c m ] σ = 50 km/s σ = 100 km/s σ = 150 km/s σ = 200 km/s b [day]10 ρ D Mm a x [ G e V / c m ] β = π /4 β = π /4v w =220 km/sv w =100 km/s FIG. 5. Upper bounds on the DM density due to DM dynamical friction as a function of the orbital period P b of a binary pulsarwith excess secular change | ˙ P xs b | (cid:46) × − and for different choices of the DM velocity dispersion σ . We considered β = π/ m = 1 . M (cid:12) and m = 0 . M (cid:12) . Right panel: m = 2 M (cid:12) and m = 5 M (cid:12) . The horizontal line denotes the value ρ DM ≈ × GeV / cm ,corresponding to a distance D ≈ as small as ˙ P int b ∼ − , the constraints shown in Fig. 5are not competitive compared to the local value of DMdensity, ρ DM ≈ . / cm . Nonetheless, Fig. 5 showsthat detecting a binary pulsar at D (cid:46)
10 pc would givestrong constraints on theoretical DM halo profiles and onthe evolution of the central BH, e.g. on the formation ofDM spikes [35].Note that the most relativistic binary pulsars have P b = O (0 .
1) day, including J1738+0333 which is the sys-tem with the most precise measurement of ˙ P xs b consistentto zero known to date [8]. As previously discussed, ouranalysis does not apply to these systems, because Eq. (4)is not satisfied. On the other hand, at least in the large- P b regime, the effect of DM dynamical friction grows lin-early with P b so the most stringent upper bounds comefrom binary systems at large orbital distance. A. Upper limit on DM density derived from thetiming of J1713+0747
Binary pulsar J1713+0747 [39] is particularly wellsuited for our analysis, due to its long orbital period, P b ≈ . | ˙ P int b | (cid:46) × − [40]. Since gravitational-wave dissipa-tion is negligible for this system ( ˙ P GW b ≈ − × − , aspredicted by the quadrupole formula [32]), ˙ P xs b ≈ ˙ P int b .For this system, dynamical friction can be larger thangravitational-wave dissipation by several orders of mag-nitude and would be the dominant contribution to theintrinsic orbital-period decay. J1713+0747 is in a nearlycircular orbit ( e ≈ × − , so our assumption of cir-cular motion at Keplerian order is accurate) and it isformed by a neutron star with m ≈ . M (cid:12) and by a white dwarf with m ≈ . M (cid:12) . Its measured inclina-tion and longitude of the pericenter are ι ≈ . ◦ and ω ≈ . ◦ , respectively [40], whereas the distance fromthe Galactic center is about D ≈ . σ ≈
150 km / s and v w ∼ (240 ±
31) km / s. Thelatter range is estimated by adding the measured trans-verse velocity of J1713+0747, v T ∼ .
28 km / s [37, 38],to the Galactic rotational velocity at the Sun location, v (cid:12) ≈
240 km / s [41].By inserting these values into Eqs. (18), (20) and (21),we can evaluate the total secular change of the orbitalperiod and estimate the maximum value of ρ DM such that |(cid:104) ˙ P DF b + ˙ P cm b + ˙ P ιb (cid:105)| (cid:46) × − , which is the value of ˙ P int b observed for J1713+0747 [40]. This condition implies theupper bound on ρ DM shown in Fig. 6. From these resultswe obtain ρ DM (cid:46) (1 − × GeV / cm , (28)where the boundaries of the range correspond to themaximum and minimum values of ρ max DM in the range α ∈ [0 , π/ β ∈ [0 , π/
2] and v w ∈ [209 , / s. Asclear from Fig. 6, the upper bound is only mildly depen-dent on the angle α and displays only small variations In principle, the values of v w , α , and β can be extracted froma measurement of the three-dimensional velocity of the binary.However, for most binary pulsars, an estimate of the radial ve-locity is not available. An exception is J1738+0333 whose three-dimensional velocity has been recently estimated [8]. For ourheuristic purposes, the estimate v w = v (cid:12) ± v T is sufficientlyaccurate. We recall that ˙ P DF b is independent of α , so the mild dependenceshown in Fig. 6 comes entirely from the two other terms, ˙ P cm b and ˙ P ιb , and mainly from the latter. in the entire range β ∈ [0 , π/
2] and v w ∈ [209 , / s.Thus, the upper bound on ρ DM is quite solid against theuncertainties on the three-dimensional velocity of this bi-nary pulsar. Finally, the upper bound (28) scales linearlywith ˙ P xs b so any future improvement on the measurementof ˙ P xs b for J1713+0747 would provide a more stringentconstraint on ρ DM . ρ D Mm a x [ G e V / c m ] v w = 209 km/sv w = 240 km/sv w = 271 km/s β [degree]02468 ρ D Mm a x [ G e V / c m ] α = π /2 α =0 FIG. 6. Upper bounds on the DM density due to DM dynami-cal friction as a function of the angle β of the DM wind derivedfrom the timing of binary pulsar J1713+0747. We adoptedthe observed values P b ≈ . e ≈ m ≈ . M (cid:12) , m ≈ . M (cid:12) , ι ≈ . ◦ , ω ≈ . ◦ and considered v w ∼ (240 ±
31) km / s. The upper (lower) panel shows thecase with α = π/ α = 0). The behavior is monotonic with α so these two cases bracket the minimum and the maximumupper bound. IV. DISCUSSION
The very fact that pulsar-timing techniques havereached the precision to be sensitive to the effect of DM inbinary pulsars is certainly intriguing. Our results showthat these effects are observable for various binaries inDM-rich environments. DM dynamical friction intro-duces a peculiar annual modulation of the orbit, similarin spirit to the one expected for the scattering rate inDM direct-detection experiments on Earth [12, 13].Current observations can be used to derive mild con-straints on the local DM density, but are limited by thefact that none of the observed binary pulsars are locatedvery close to the Galactic center, where the DM densityis expected to peak. The next Square Kilometer Array(SKA) telescope [42] can improve the precision of cur-rent measurements by 2 orders of magnitude, potentiallyprobing DM densities lower than the DM density in theSolar System. SKA is also expected to increase the num-ber of observed pulsars by a factor of 10, thus enhanc-ing the chances of observing a binary pulsar in regions of high DM density. Our analysis suggests that precisemodels of the effects of DM on the orbital motion wouldbe mandatory to reach high accuracy in timing modelsfor such binaries.In this paper we presented an exploratory study thatis valid only for binaries with long orbital period. Formore relativistic binaries with P b (cid:46) . P b in some systems can be used to discrimi-nate the effect of DM from other mechanisms [23, 29–31]that might cause comparable changes in the intrinsic or-bital period. A more sophisticated analysis can improvethe constraints derived in this paper and would be crucialto devise detection tools, rather than to simply put upperbounds on the DM density. We hope that the prospectsraised by our analysis will encourage further research inthis direction. ACKNOWLEDGMENTS
We are indebted to Norbert Wex and Dick Manchesterfor very instructive correspondence, to Matteo Cadeddu,Vitor Cardoso and Leonardo Gualtieri for interesting dis-cussions, and to two anonymous referees for several valu-able comments. This work was supported by the Eu-ropean Community through the Intra-European MarieCurie Contract No. AstroGRAphy-2013-623439, byFCT-Portugal through Project No. IF/00293/2013, bythe NRHEP 295189 FP7-PEOPLE-2011-IRSES Grant,and by the COST Action MP1304 “NewCompStar”.
Appendix A: On the dynamical friction formula
Chandrasekhar’s dynamical friction formula for a sin-gle object moving through a homogeneous density dis-tribution [14] is derived under a number of assumptions.Here we repeat the standard derivation and discuss thevalidity of certain assumptions.Let us consider a star with mass m and velocity v moving through a homogeneous density distributionmade of light particles with mass m DM (cid:28) m . Due tothe encounter with a particle of mass m DM moving withvelocity u , the velocity of the star is incremented by ∆ v (cid:107) and ∆ v ⊥ , in the directions parallel and perpendicular to v , respectively. The total incremental velocity sufferedby the star can be obtained by summing all contributionsof the particles with mass m DM over an interval of timewhich is much longer than the typical interaction timebut also much shorter than the time scale over which v changes. As expected by symmetry arguments, whensummed over a large number of encounters (cid:80) ∆ v ⊥ van-ishes, while the incremental change in the direction ofmotion reads [15] d v dt = − π ( m + m DM ) m DM G v v (cid:90) ∞ d uf ( u ) Q ( u ) , (A1)where m DM f ( u ) is the density distribution and Q ( u ) = log (cid:2) (1 + q ( v + u ) )(1 + q ( v − u ) ) (cid:3) u < v log (cid:2) (1 + 16 q u (cid:3) − u = v log q ( u + v ) q ( u − v ) − vu u > v , with q = b max / [ G ( m + m DM )] and b max being the char-acteristic size of the medium. Assuming q ( v ± u ) (cid:29) d v dt = − π G ( m + m DM ) m DM v v × (cid:20)(cid:90) v d uf ( u ) (cid:18) log( qv ) + log v − u v (cid:19) + (cid:90) ∞ v d uf ( u ) (cid:18) − vu + log u + vu − v (cid:19)(cid:21) . (A2)If the velocity dispersion σ is large, σ (cid:29) v , the aboveequation reduces to d v dt ∼ − π G ( m + m DM ) m DM v v × (cid:20) log( qv ) (cid:90) v d uf ( u ) + 23 v (cid:90) ∞ v d u f ( u ) u (cid:21) . (A3)Equation (A3) contains two contributions to the vis-cous drag on m , namely a term arising from slow par-ticles with velocity smaller than v and a term com-ing from fast particles with velocity larger than v . Asnoted in Ref. [24], these two terms are roughly of thesame order. For an isotropic, Maxwellian distribution m DM f ( u ) = ρ DM (2 πσ ) / e − u / (2 σ ) , these two contributionscombine to give d v dt ∼ − √ π λ ρ DM G ( m + m DM ) σ v . (A4)where λ = log( qσ ) is the Coulomb logarithm. Equa-tion (A4) coincides with the σ (cid:29) v limit of Chan- drasekhar’s formula [14, 15] d v dt = − πλ ρ DM G ( m + m DM ) v (cid:18) erf ( x ) − x √ π e − x (cid:19) v , (A5)where x = v/ ( √ σ ). Although it is usually consideredthat dynamical friction on a single body is due to am-bient particles moving slower than the perturber, in factEqs. (A2) and (A3) show that faster particles also con-tribute to the final result.In a binary system, the situation is actually opposite tothe case of dynamical friction on a single perturber [24].In that case the contribution from faster particles is dom-inant because such encounters have an impact parame-ter smaller than the orbital distance and interact witheach component of the binary system separately. On theother hand, slower particles interact with the binary asa whole through tidal forces, only perturbing its centerof mass [24, 25]. Interestingly, for a Maxwellian distri-bution the two effects are equivalent but multiplied bytwo different Coulomb logarithms, one (arising from anintegral over the impact parameter) is associated to slowencounters, whereas the other (arising from an integralover velocity space) is associated to fast encounters. Inthis case the final result is still described by Eq. (A5)with an effective Coulomb logarithm λ . This property isa great advantage for our analysis and is valid only forMaxwellian or nearly Maxwellian distributions, as thoseassumed in the main text. As noted in the main text ourresults, based on a superposition of Eq. (A5) for each sin-gle object of the binary system, are in exact agreementwith those derived by Gould [24] in the large- σ limit andwithin a completely different framework. This agreementgives further support to the validity of our analysis.Nonetheless, the final formula (A5) is obtained fromEq. (A1) after various approximations [15] whose validityis unclear (for example, while q ( v + u ) (cid:29) q ( v − u ) is a small quantity when u ∼ v ).Thus, it is relevant to compare the exact numerical resultagainst the approximate one typically used. We makethis comparison in Fig. 7. The analytical formula (A5) -2 -1 v/ σ nu m e r i c a l / ana l y t i c a l λ ~16.1 λ ~13.8 FIG. 7. Ratio between Eq. (A1) and Eq. (A5) as a functionof v . The agreement of the exact result with the analyticalformula (A5) is very good for any v < σ . v < σ and deviates at most by ∼
35% when v ∼ σ .In the main text we focus on the regime (4), where theanalytical formula is in very good agreement with theexact results. Even for mildly relativistic pulsar binarieswith v ∼ σ ∼
300 km / s, Chandrasekhar’s formula isaccurate roughly within 6%. Appendix B: Computation of ˙ P cm b Our starting points are Eq. (7) and Eq. (19). To firstorder, Eq. (7) can be written as ˙ V = a η v + a v w . We de-compose the DM wind vector as v w = v w cos α cos β e x + v w sin α sin β e y + v w cos β e z , where e x , e y and e z are the unit vectors in the orbital frame. The latter are relatedto the unit vectors e X , e Y , e Z in the fundamental frameby a sequence of rotations (cf. Ref. [22] for details). Wefurther decompose v = v n n + v λ λ + v z e z . To zeroth orderand in the circular case, v = v λ .The final expression Eq. (20) can be obtained by mak-ing use of the relations λ · e Z = sin ι cos(Ω t + ω ) , (B1) e x · e Z = sin ι sin ω , (B2) e y · e Z = sin ι cos ω , (B3) e z · e Z = cos ι , (B4)which are derived in Sec. 3.5.2 of Ref. [22]. [1] G. Bertone, ed., Particle Dark Matter (Cambridge Uni-versity Press, 2010).[2] R. A. Flores and J. R. Primack, Ap.J. Letters , L1(1994), astro-ph/9402004.[3] B. Moore, Nature , 629 (1994).[4] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, andJ. Hessels, Nature , 1081 (2010), arXiv:1010.5788[astro-ph.HE].[5] J. Antoniadis, P. C. C. Freire, N. Wex, T. M. Tauris, R. S.Lynch, M. H. van Kerkwijk, M. Kramer, C. Bassa, V. S.Dhillon, T. Driebe, J. W. T. Hessels, V. M. Kaspi, V. I.Kondratiev, N. Langer, T. R. Marsh, M. A. McLaugh-lin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J. vanLeeuwen, J. P. W. Verbiest, and D. G. Whelan, Science (2013), 10.1126/science.1233232.[6] R. Hulse and J. Taylor, Astrophys.J. , L51 (1975).[7] A. Lyne, M. Burgay, M. Kramer, A. Possenti, R. Manch-ester, et al. , Science , 1153 (2004), arXiv:astro-ph/0401086 [astro-ph].[8] P. C. C. Freire, N. Wex, G. Esposito-Far`ese, J. P. W. Ver-biest, M. Bailes, B. A. Jacoby, M. Kramer, I. H. Stairs,J. Antoniadis, and G. H. Janssen, MNRAS , 3328(2012), arXiv:1205.1450 [astro-ph.GA].[9] M. Kramer, 2, D. C. Backer, J. M. Cordes, T. J. W.Lazio, B. W. Stappers, and S. Johnston,
Proceed-ings, IAU Symposium on From Planets to Dark Energy:The Modern Radio Universe (MRU 2007) , New Astron.Rev. , 993 (2004), [PoSMRU,020(2007)], arXiv:astro-ph/0409379 [astro-ph].[10] N. Yunes and X. Siemens, Living Rev.Rel. , 9 (2013),arXiv:1304.3473 [gr-qc].[11] E. Berti et al. , Class. Quant. Grav. , 243001 (2015),arXiv:1501.07274 [gr-qc].[12] A. K. Drukier, K. Freese, and D. N. Spergel, Phys. Rev.D , 3495 (1986).[13] K. Freese, M. Lisanti, and C. Savage, ArXiv e-prints(2012), arXiv:1209.3339 [astro-ph.CO].[14] S. Chandrasekhar, Astrophys.J. , 255 (1943).[15] S. Chandrasekhar, Rev. Mod. Phys. , 383 (1949).[16] J. Binney and S. Tremaine, Galactic Dynamics (Prince-ton University Press, Princeton, NJ, 1987).[17] S. W. Randall, M. Markevitch, D. Clowe, A. H. Gon-zalez, and M. Bradac, Astrophys. J. , 1173 (2008), arXiv:0704.0261 [astro-ph].[18] D. N. C. Lin and S. Tremaine, Astrophysical Journal , 364 (1983).[19] H. Kim and W.-T. Kim, Astrophys. J. , 432 (2007),arXiv:0705.0084 [astro-ph].[20] A. Just, F. M. Khan, P. Berczik, A. Ernst, andR. Spurzem, MNRAS , 653 (2011), arXiv:1009.2455[astro-ph.CO].[21] J. D. Bekenstein and R. Zamir, Astrophysical Journal , 427 (1990).[22] E. Poisson and C. Will,
Gravity: Newtonian, Post-Newtonian, Relativistic (2014).[23] D. R. Lorimer and M. Kramer,