Dark-matter dynamical friction versus gravitational-wave emission in the evolution of compact-star binaries

TL;DR

This paper investigates the impact of dark matter dynamical friction on the orbital evolution of compact-star binaries, comparing it with gravitational wave emission, and explores how this effect varies with Galactic position and dark matter profiles.

Contribution

It provides the first detailed analysis of dark matter dynamical friction effects on binary evolution, considering multiple DM profiles and their influence on orbital decay or widening.

Findings

GW emission dominates over DMDF at typical Galactic positions.

DMDF can become significant at larger orbital periods or closer to the Galactic center.

Potential to constrain dark matter properties through binary orbital observations.

Abstract

The measured orbital period decay of compact-star binaries, with characteristic orbital periods $\sim 0.1$~days, is explained with very high precision by the gravitational wave (GW) emission of an inspiraling binary in vacuum. However, the binary gravitational binding energy is also affected by an usually neglected phenomenon, namely the dark matter dynamical friction (DMDF) produced by the interaction of the binary components with their respective DM gravitational wakes. The entity of this effect depends on the orbital period and on the local value of the DM density, hence on the position of the binary in the Galaxy. We evaluate the DMDF produced by three different DM profiles: the Navarro-Frenk-White (NFW), the non-singular-isothermal-sphere (NSIS) and the Ruffini-Arg\"uelles-Rueda (RAR) profile based on self-gravitating keV fermions. We first show that indeed, due to their Galactic position, the GW emission dominates over the DMDF in the NS-NS, NS-WD and WD-WD binaries for which measurements of the orbital decay exist. Then, we evaluate the conditions under which the effect of DMDF on the binary evolution becomes comparable to, or overcomes, the one of the GW emission. We find that, for instance for $1.3$--$0.2$ $M_\odot$ NS-WD, $1.3$--$1.3$~$M_\odot$ NS-NS, and $0.25$--$0.50$~$M_\odot$ WD-WD, located at 0.1~kpc, this occurs at orbital periods around 20--30 days in a NFW profile while, in a RAR profile, it occurs at about 100 days. For closer distances to the Galactic center, the DMDF effect increases and the above critical orbital periods become interestingly shorter. Finally, we also analyze the system parameters for which DMDF leads to an orbital widening instead of orbital decay. All the above imply that a direct/indirect observational verification of this effect in compact-star binaries might put strong constraints on the nature of DM and its Galactic distribution.

Tables

Table 1. Table 1: Intrinsic orbital decays for several binary systems in the Galaxy as well as the ones predicted by GR and DM dynamical friction. There, it is also shown the values of mass binaries, orbital periods and distances measured from the Galactic center. This information is taken completely from Table 1 . in Ref.  Antoniadis ( 2015 ) and references therein. For updated values of masses of neutron stars see Özel and Freire ( 2016 ) . We have simply added the last row for the WD-WD binary and the last two columns to show the orbital decay predicted by DMDF for the NFW profile and the RAR model.

Name	Type	$m_p$ [$M_{\odot }$]	$m_c$ [$M_{\odot }$]	$P_b$ [days]	d [kpc]	${\dot{P}}_b^{\mathrm{int}}$ [${10}^{-12}$]	${\dot{P}}_b^{\mathrm{GW}}$ [${10}^{-12}$]	${\dot{P}}_{b,NFW}^{\mathrm{DF}}$ [${10}^{-21}$]	${\dot{P}}_{b,\mathrm{RAR}}^{\mathrm{DF}}$ [${10}^{-21}$]
J0737-3039	NS-NS	1.3381(7)	1.2489(7)	0.104	1.15(22)	-1.252(17)	-1.24787(13)	-10.498	-7.860
B1534+12	NS-NS	1.3330(4)	1.3455(4)	0.421	0.7	-0.19244(5)	-0.1366(3)	-244.166	-27.827
J1756-2251	NS-NS	1.312(17)	1.258(17)	0.321	2.5	-0.21(3)	-0.22(1)	-0.271	-20.695
J1906+0746	NS-NS	1.323(11)	1.290(11)	0.166	5.4	-0.565(6)	-0.52(2)	-2.655	-11.176
B1913+16	NS-NS	1.4398(2)	1.3886(2)	0.325	9.9	-2.396(5)	-2.402531(14)	-7.942	-17.747
B2127+11C999This binary is located in the globular cluster M15 Anderson et al. (1990). However we have made here a simple estimation of the DMDF effect assuming that the DM local density in its location does not change abruptly within the globular cluster, which may not be the case. This point is better discussed in footnote 10. For a comprehensive list of all known binaries in globular clusters see http://www.naic.edu/ pfreire/GCpsr.html and references therein.	NS-NS	1.358(10)	1.354(10)	0.333	10.3(4)	-3.961(2)	-3.95(13)	-8.083	-17.0154
J0348+0432	NS-WD	2.01(4)	0.172(3)	0.104	2.1(2)	-0.273(45)	-0.258(11)	-0.399	-1.514
J0751+1807	NS-WD	1.26(14)	0.13(2)	0.263	2.0	-0.031(14)	—	-1.022	-2.587
J1012+5307	NS-WD	1.64(22)	0.16(2)	0.60	0.836(80)	-0.15(15)	-0.11(2)	-3.404	-7.343
J1141-6545	NS-WD	1.27(1)	1.02(1)	0.20	3.7	-0.401(25)	-0.403(25)	-3.578	-11.469
J1738+0333	NS-WD	1.46(6)	0.181(7)	0.354	1.47(10)	-0.0259(32)	-0.028(2)	-2.120	-4.379
WDJ0651+2844	WD-WD	0.26(4)	0.50(4)	0.008	1	-9.8(28)	-8.2(17)	-0.014	-0.207

Table 2. Table 2: This table displays theoretical predictions of orbital periods in days at which P ˙ b DF superscript subscript ˙ 𝑃 𝑏 DF \dot{P}_{b}^{\rm DF} , computed by the indicated DM density profiles, equates P ˙ b GW superscript subscript ˙ 𝑃 𝑏 GW \dot{P}_{b}^{\rm GW} predicted by general relativity for different binary systems.

DM Profile	NS-WD	NS-NS	WD-WD
NFW	18	30	25
RAR	120	130	150
NSIS	80	90	70

Equations

$${\dot{P}}_b^{\mathrm{GW}}=-\frac{192\pi }{5}{\left(\frac{2\pi G\mathcal{M}}{c^3{P}_b}\right)}^{5/3},$$

$${\text{f}}_{fr,i}$$

$$+{\int }_{{\tilde{v}}_i}^{v_{\mathrm{esc}}}{d}^{3}uf(u)\left[\ln \left(\frac{u+{\tilde{v}}_i}{u-{\tilde{v}}_i}\right)-2\frac{{\tilde{v}}_i}{u}\right])\frac{{\tilde{\text{v}}}_i}{{\tilde{v}}_i^3},$$

$$\log \Lambda \equiv \log \left(\frac{b_{\max }}{b_{\min }}\right).$$

$$b_{\max }\approx a,b_{\min }=\max (r_h,R_A),$$

$$R_{A,i}=\frac{G{m}_i}{{\tilde{v}}_i^2},$$

$$f(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}\left[{\int }_0^{\mathcal{E}}\frac{d^2\rho }{d{\Psi }^2}\frac{d\Psi }{\sqrt{\mathcal{E}-\Psi }}+\frac{1}{\sqrt{\mathcal{E}}}{\left(\frac{d\rho }{d\Psi }\right)}_{\Psi =0}\right],$$

$$f^{\mathrm{MB}}(u)=\frac{n_0}{(2\pi {\sigma }^2{)}^{3/2}}\exp \left(-\frac{u^2}{2{\sigma }^2}\right),$$

$$f_{c}(p)=\frac{gm^{3}}{h^{3}}\left\{\begin{array}[]{ll}\frac{1-e^{(\epsilon-\epsilon_{c})/kT}}{e^{(\epsilon-\mu)/kT}+1}&\epsilon\leqslant\epsilon_{c},\\ 0&\epsilon>\epsilon_{c}.\end{array}\right.$$

$$\varphi (r)=4\pi G\left[\frac{1}{r}{\int }_0^{r}d{r}^{\prime }{r}^{\prime 2}\rho (r^{\prime })+{\int }_r^{\infty }d{r}^{\prime }{r}^{\prime }\rho (r^{\prime })\right].$$

$$\rho (r)=\frac{{\rho }_c}{(r/r_s)(1+r/r_s{)}^2},$$

$$\rho (r)=\frac{{\rho }_0}{1+(r/r_0{)}^2},$$

$$\rho$$

$$P$$

$$\frac{1}{\rho (r)}\frac{d(\rho (r){\sigma }_r^2)}{dr}+2\frac{\beta {\sigma }_r^2}{r}=-\frac{d\varphi (r)}{dr}=-V_c^2,$$

$$v_c^2=-{\sigma }_r^2\left(\frac{d\ln \rho }{d\ln r}+\frac{d\ln {\sigma }_r^2}{d\ln r}+2\beta \right),$$

$$\dot{\text{v}}=-\frac{GM}{r^3}\text{r}+\text{f},$$

$$\dot{a}=2\sqrt{\frac{r_0^3}{GM}}S(t),$$

$$\dot{e}=2\sqrt{\frac{r_0}{GM}}[R(t)\sin ({\Omega }_{0}t)+2S(t)\cos ({\Omega }_{0}t)],$$

$$\dot{i}=2\sqrt{\frac{r_0}{GM}}W(t)\cos ({\Omega }_{0}t+\omega ),$$

$$\dot{\Omega }=\frac{1}{\sin i}\sqrt{\frac{r_0}{GM}}W(t)\sin ({\Omega }_{0}t+\omega ),$$

$$S(t)=a_1\eta v-a_2{v}_w\sin \beta \sin ({\Omega }_{0}t-\alpha ),$$

$$R(t)=a_2{v}_w\sin \beta \cos ({\Omega }_{0}t-\alpha ),$$

$$W(t)=a_2{v}_w\cos \beta .$$

$$\frac{\dot{P_b}}{P_b}=\frac{3}{2}\frac{\dot{a}}{r_0}.$$

$$\dot{P_b}(t)=3P_b[a_1\eta -a_2\Gamma \sin \beta \sin ({\Omega }_{0}t-\alpha )].$$

$$⟨{\dot{P}}_b⟩=\frac{1}{P_b}{\int }_0^{P_b}\dot{P_b}(t)dt.$$

$$b_i=\frac{1}{\rho (r)\log {\Lambda }_i}\frac{I_i}{{\tilde{v}}_i^3},$$

$$a_1=-(A_1{b}_1+A_2{b}_2),a_2=\frac{1}{2}(A_1{b}_1{\Delta }_{+}+A_2{b}_2{\Delta }_{-}),$$

$$\frac{d^2\bar{\rho }}{d{\bar{\Psi }}^2}=\frac{d}{d\bar{r}}\left(\frac{\bar{\rho }\prime (\bar{r})}{\bar{\Psi }\prime (\bar{r})}\right)\frac{d\bar{r}}{d\bar{\Psi }}.$$

$$f(\mathcal{E})=\frac{{\bar{\rho }}_0}{(2\pi {\sigma }^2{)}^{3/2}}\exp [\mathcal{E}/{\sigma }^2],$$

Full text

DM: density profiles, velocity distribution function, halo- Galaxy: dynamical friction, Pulsars: Binaries, orbital period decay

I Introduction

Compact-star binaries composed of neutron stars (NSs) and/or white dwarfs (WDs) have turned out to be rich laboratories of physics and astrophysics that allow to test fundamental theoretical predictions. In particular, NS-NS binaries have served to prove the existence of gravitational waves (GWs) Hulse and Taylor (1975) and the motion of matter and photons in the strong gravitational fields Kramer et al. (2006), as well as other phenomena Lorimer and Kramer (2012). These latter aspects are of special interest in tests of general relativity and alternative theories of gravity Kramer et al. (2006); Antoniadis et al. (2013).

The orbital motion of such systems also offers the possibility of analyzing further effects. An interesting physical situation arises when the orbiting object moves through an extended medium which is formed, for instance, from the mass loss of the binary companion. This interaction can be thought as a drag force exerted by the circumbinary medium on the object in question, perturbing thereby its Keplerian orbital motion Alexander et al. (1976). This dynamical friction produced by the gravitational drag-force has been also studied in the context of different astrophysical phenomena such as mergers of star clusters, galaxies, and even galaxy clusters, to the inspiral of dwarf galaxies within dark-matter halos and the orbital evolution of massive black hole (BH) binaries in a stellar medium Binney and Tremaine (2008). Thus, dynamical friction plays an important role in the orbital evolution of many astrophysical systems. In a pioneering work, S. Chandrasekhar Chandrasekhar (1943) calculated the dynamical friction force on a massive object traversing an infinite homogeneous collisionless background (representing the surrounding star neighbors).

It is thus natural to expect that a binary system moving through the galaxy can also experience a dynamical friction caused by collisionless DM particles, namely DM dynamical friction (hereafter DMDF), particularly in DM-dominated regions, as at the outer part of the Galactic halo and near the Galactic center Strigari (2013). The perturbed orbital motion may lead thus to interesting observable effects in the secular evolution of the orbital period. An interesting proposal was advanced in Ref. Pani (2015) on the possibility of inferring constraints to the DM density by determining the above DM effect on the orbital motion of binaries (see also the pioneering work by Bekenstein & Zamir Bekenstein and Zamir (1990), for a general discussion of collisionless background types as well as in the context of DM). They showed that the change in the orbital period could be due to the dynamical friction force exerted by the DM background on the binary. In that work, this effect was used to put an upper bound on the DM density in a given location of the Galaxy, independently of the density profile or the nature of the DM particles. It can be shown, however, that this upper limit is indeed fulfilled by any DM density profile consistent with the outer halo properties of the Milky Way. Thus, we explore in this work the dependence of the orbital period decay by DMDF on the different binary parameters and also on the DM density profile, in order to identify all possible physical situations suitable for an observational verification of the DMDF effect. For doing this we obtain DM profiles fulfilling definite Galactic-halo observables such as the escape velocity, the velocity dispersion and the one-halo scale length parameters. The velocity distribution function and the DM density profile are, as we shall show below, crucial elements in the dynamical friction force estimation.

It is known that the DM in the outer part of our Galaxy is well described by a classical Maxwell-Boltzmann distribution, e.g. by a non-singular isothermal (hereafter NSIS) profile Binney and Tremaine (2008). However, depending on the DM nature (e.g. particle type), the DM density distribution can deviate from the classical Maxwell-Boltzmann behavior towards the inner regions of the Galaxy. This implies that the DMDF effect will depend according to the phase-space density consistent with the DM particle nature. We shall consider, for the sake of comparison, three DM models: 1) the NSIS profile, 2) the Navarro-Frenk-White (NFW) profile Navarro et al. (1996), and 3) the recently introduced Ruffini-Argüelles-Rueda (RAR) model Ruffini et al. (2015); Argüelles et al. (2016).

The RAR model is based on a self-gravitating system of massive (keV) fermions in thermodynamic equilibrium. The density profile of the RAR model exhibits a core-halo structure which allows to explain the DM distribution in galactic halos from dwarfs to big spirals, and predicts at the same time the presence of a DM high density core Ruffini et al. (2015). Under this approach and following the more realistic distribution function including violent relaxation processes (Chavanis, 2004) and the escape velocity of particles, the Fermi-Dirac distribution function was subsequently introduced to describe the finite size of halos. In the case of the Milky Way, such a DM core can explain the observed dynamics near the Galactic center Sgr A* without invoking a central supermassive black hole for fermion masses in the range 48 keV $\lesssim mc^{2}\lesssim 345$ keV Argüelles et al. (2016) 111See also Ref. Gómez et al. (2016) for the gravitational lensing properties of the RAR profile..

Having established the DM density profiles we shall analyze, we now describe the structure of this work. We start by discussing in section II the effects which are commonly assumed to produce a change of the orbital period of binaries, putting special attention evidently to the one produced by GW emission. We analyze in section III the dynamical friction force and its main ingredients for the case when it is produced by DM and when it acts on binary systems. We analyze in section IV the perturbation effect of DMDF on the orbital motion of the pulsar and reproduce some general results presented in Pani (2015). Furthermore, we introduce Galactic-halo observables in order to generalize the prescription presented in Pani (2015) and present thus a more realistic estimation of dynamical friction effects. Finally, we present in section V the numerical results of $\dot{P}_{b}$ as a function of the radial position, the DM wind velocity and the orbital period. This latter computation leads us to compare directly the $\dot{P}_{b}$ due to GW emission to that given by DMDF. In section V we summarize our results and present a general discussion.

II Binary systems and orbital period decay by gravitational waves

The precise pulsar timing measurements allow us to detect, with a high accuracy, tiny orbital effects which thus require a precise theoretical description of the orbital motion Hulse and Taylor (1975). In the weak field regime (Newtonian approach), the binary motion of pulsar is simply described by the Kepler laws. However, relativistic and strong-field effects in the orbital motion should be taken into account in the vicinity of a close-orbit binary pulsar Kramer et al. (2006). These relativistic effects can be described, for the known binaries, with sufficient accuracy in terms of the called post-Keplerian parameters that account for departures from Newtonian Keplerian dynamics owing e.g. to the GW emission, time delay caused by the curvature of space-time near the pulsar (Shapiro delay), and relativistic time dilation Poisson and Will (2014). There exists a variety of effects that affect the orbital period stability and they can be, roughly speaking, classified in two large groups: kinematic and intrinsic to the system. The former include the effects of a secular increase due to the Galactic gravitational potential, secular acceleration resulting from the pulsar’s transverse velocity (proper motion of the pulsar) and the cluster’s gravitational field; while the latter is related to “local” effects in the system as mass loss either from the pulsar or its companion and the GW emission among others222For a more detail description of possible effects on the observed period decay see Refs. Damour and Taylor (1991); Lorimer and Kramer (2012).. After subtracting kinematic effects from the observed change of the orbital period, the remaining intrinsic period decay has been shown to be explained by the GW emission predicted by general relativity of an inspiraling binary in vacuum.

The orbital period decay owing to the GW emission of a binary spiraling in circular orbits is given by

[TABLE]

where $G$ is the gravitational constant, $P_{b}$ is the orbital period, $\mathcal{M}=(m_{p}m_{c})^{3/5}M^{-1/5}$ is the so-called chirp mass and $M=m_{p}+m_{c}$ is the total mass, with the subscripts $p$ and $c$ denoting the primary component and its companion, respectively.

The theoretical prediction of general relativity given by Eq. (1) was first verified with the observed intrinsic orbital period decay of the famous Hulse-Taylor binary pulsar PSR B1913+16 Hulse and Taylor (1975), which is explained with an accuracy of 99,8%. Later on, additional successful verifications in other relativistic NS-NS and NS-WD binaries have been made and with even higher accuracy. We refer the reader to Ref. Antoniadis (2015) for a review on this subject and also Table 1.

As we have mentioned, the above orbital period decay by GW emission is calculated under the assumption of binary motion in empty space. We shall explore below the effect of the presence of DM background on the orbital motion via dynamical friction, i.e. by DM gravitational drag. We shall infer the predicted orbital period time derivative by this phenomenon to then compare it with the one produced by the GW emission.

III Dynamical friction force and its main ingredients

Dynamical friction has been widely used to account for the drag force when an object is moving through a collisionless medium of field particles. This drag induces a wake of medium particles on the object with a characteristic overdensity proportional to its mass Binney and Tremaine (2008). In his seminal work, Chandrasekhar Chandrasekhar (1943) computed the dynamical friction force onto an object that move in an infinite homogeneous stellar medium obeying a Maxwell-Boltzmann velocity distribution, taking into account only the contribution of the field particle velocities smaller with respect to the object’s velocity. However, the dynamical evolution of many astrophysical systems is driven by dynamical friction in a more realistic way Binney and Tremaine (2008). We consider the drag force, $\textbf{f}_{fr,i}$, experienced by a test body of mass $m_{i}\gg m$, being $m$ the DM particle mass, and with orbital velocity $v_{i}$ moving through the DM background with velocity distribution function $f(u)$ Chandrasekhar (1943); Binney and Tremaine (2008):

[TABLE]

where the integral in the first term, accounts for low velocity contributions (fraction of particles moving slower than the object) while the integral in the second term, refer to the faster particles, limited by the escape velocity $v_{\rm esc}$ according to the Galactic gravitational potential.333It has been recently shown that the incorporation of the tidal radius into the background system can produce interesting features in infalling satellites in large cored galaxies Petts et al. (2016). $b_{\rm max}$ is the maximum impact parameter defined below in Eq. (4). The above equation takes into account the orbital velocity of each object with respect to the DM wind relative to the center of mass of the binary system: $\tilde{\textbf{v}}_{i}=\textbf{v}_{i}+\textbf{v}_{w}$, with $\textbf{v}_{w}=v_{w}(\cos{\alpha}\sin{\beta},\sin{\alpha}\sin{\beta},\cos{\beta})$ and $\beta$ and $\alpha$ being the angles between the wind velocity vector and the perpendicular axis of the binary orbital plane and the projection of the wind velocity vector with an axes lying in the orbital plane, respectively. There are at least two different cases of wind velocities: bound and unbound binaries to the galaxy potential. In the former the DM wind velocity can be assumed as the negative of the binary circular velocity with respect to the galactic center $\textbf{v}_{w}=-\textbf{v}_{\rm rot}$. The latter case occurs often in binaries with NS components in which the system received a high kick velocity from the supernova event Lyne and Lorimer (1994). For high kick velocities the binary circular velocity with respect to the galactic center can be neglected Postnov and Yungelson (2006) and we can assume $\textbf{v}_{w}=-\textbf{v}_{T}$, where $\textbf{v}_{T}$ is the transversal velocity of the system. For intermediate kicks, the system can remain bound and we can consider, in a more general case, $\textbf{v}_{w}=\textbf{v}_{\rm rot}+\textbf{v}_{T}$. Thus, we shall consider the value of $v_{w}$ as a free parameter that can assume values ranging from 10 km s*-1* all the way to 1000 km s*-1* following the above discussion. There is the additional possibility for the binary components to experience an intrinsic DM wind. However, up to the best of our knowledge, there is no observational evidence of an intrinsic rotation of the DM with respect to the Galactic center and thus we do not consider it in our estimates.

It is important also to mention that the condition $L/a\ll 1$, where $L$ is the size of the component’s wake and $a$ the orbital separation, must be fulfilled in order that Eq. (2) becomes linearly applicable to each binary component Bekenstein and Zamir (1990); Binney and Tremaine (2008). Since $L$ is of the order of the radius of the sphere of gravitational influence of each component – see Eq. (5) below – this means that we are limited to binary systems with orbital velocities smaller than the velocity dispersion of the DM background. Namely, we deal with binary systems with sufficiently large orbital periods (small orbital compactness) so that each binary component does not interact with its respective companion’s wake. Furthermore, we treat the binary system as composed of point masses no matter their internal structure. Thus, we can apply this approach under the above conditions to binary systems such as NS-NS/NS-WD Antoniadis (2015) and WD-WD Hermes et al. (2012), or any other possible binary system of astrophysical interest.

We proceed now to introduce the most relevant ingredients entering into the computation of the dynamical friction force on the binary system. This analysis allow us to establish more accurately our system in terms of Milky Way galactic observables and to characterize more realistically the DM density properties.

III.1 The Coulomb logarithm

The Coulomb logarithm in the Chandrasekhar’s dynamical friction formula accounts for the finite size of the system and is defined as the ratio of the maximum and minimum impact parameters for encounters, respectively $b_{\rm max}$ and $b_{\rm min}$, i.e.

[TABLE]

It is assumed typically that $b_{\rm max}$ is of the order of the size of the system, and $b_{\rm min}$ is defined as the impact parameter for a 90∘ deflection Binney and Tremaine (2008)

[TABLE]

where $b_{\rm max}$ can be taken as the effective size of the system (the binary orbital separation) and $r_{h}$ is the half-mass radius of the subject system. This is the radius that contains the body’s half-mass and should be taken as $b_{\rm min}$ in the case it be an extended body. However it does not correspond to the present case. We instead adopt $b_{\rm min}=R_{A}$, where $R_{A}$ is the radius at which a particle of the surrounding medium is affected by the sphere of gravitational influence of the test body, namely:

[TABLE]

being $\tilde{v}_{i}$ the relative velocity of the object with respect to the DM wind velocity as we defined above.

We can see from here that dynamical friction force is determined by the local distribution of matter producing the wake around each object. This also establishes the characteristic size of the wake. It is here assumed that $b_{\rm max}\gg b_{\rm min}$ and $b_{\rm max}$ is set to be the length scale over which the density can be assumed to be constant for a given system at fixed radial position. It is important to note that the choices of the impact parameters are somewhat arbitrary. However, we guarantee that the condition $\Lambda\gg 1$ is satisfied.

As an example we plot in Fig. 1 the Coulomb logarithm as a function of the wind velocity for a $1.3+0.2~{}M_{\odot}$ NS-WD binary with orbital period $P_{b}=100$ days. We stress that the Coulomb logarithm does not change with $\beta$, we choose however $\beta\neq 0$ to perform, in a more general way, a study of the orbital period decay. We also note that, for $v_{w}\geq 80$ km s*-1* there are not large differences between the Coulomb logarithm for each object. However, we will take into account these small differences for accuracy even though we consider, in some cases, large values for the wind velocity.

III.2 Velocity distribution function

The evolution of a collisionless self-gravitating system is determined by the Vlasov-Poisson equation that sets the conservation of the phase-space density Binney and Tremaine (2008). This distribution function fully specifies the dynamic of a colissionless system. For instance, for spherical systems, the mass density is proportional to $\int d^{3}vf$. It is also possible to derive the distribution function of a collisionless system for a given self-consistent density profile $\rho$ following the Eddington’s formula Eddington (1916)

[TABLE]

where we have introduced the relative potential and binding energy (per unit mass) defined respectively as: $\Psi=-\Phi+\Phi_{0}$ and $\mathcal{E}=-E+\Phi_{0}=\Psi-\frac{1}{2}v^{2}$. For a spherical system with an isotropic velocity dispersion, the phase space distribution function of dark halos depends only on the energy and not on the angular momentum. The above formula is particularly useful when we seek for a distribution function to associate with a density profile obtained from other methods. We shall apply this procedure in the appendix, to the NSIS and to the phenomenological NFW profile to validate the approximation of considering, within our estimations, the Maxwell-Boltzmann distribution for these both profiles. Significant but very small differences appear between the Maxwell-Boltzmann distribution and the distribution functions associated to the NFW and NSIS profile at nearby unbound energies, as can be seen in Fig. 7. In addition, the unbound energy ($\mathcal{E}=0$) permits, in fact, to define the escape velocity $v_{esc}=\sqrt{2|\Psi|}$. We shall see that the contribution of particles moving faster than the object and limited by the escape velocity, do not contribute substantially to the dynamical friction force. This consequence supports the fact of considering the Maxwell-Boltzmann distribution to describe the velocity distribution for the aforementioned profiles. The main motivation of this approach is then, due to the numerical facilities that the simple Maxwell-Boltzmann distribution provides in the computation of the dynamical friction force.

Accordingly, for the sake of comparison, let us assume then that the virialized NFW and NSIS halos, follow the Maxwell-Boltzmann distribution function

[TABLE]

where $n_{0}$ is the particle number such that $\rho=n_{0}m$ and $\sigma$ is the velocity dispersion which is defined in terms of the DM gravitational potential through the Jeans Eq. (14).

For the RAR model, we consider self-consistently a Fermi-Dirac distribution function with energy cutoff $\epsilon_{c}$ to describe the velocity distribution of self-gravitating halos in thermodynamic equilibrium Ruffini et al. (2015)444See also Ref. Chavanis (2004) for a general discussion about the conditions under which statistical equilibrium state is reached.

[TABLE]

Here $\epsilon=\sqrt{c^{2}p^{2}+m^{2}c^{4}}-mc^{2}$ is the particle kinetic energy, $m$ is the particle mass, $\mu$ is the chemical potential (with the particle rest mass subtracted off), $T$ is the temperature, $k$ is the Boltzmann constant. The quantity $g$ denotes as usual the particle spin degeneracy ($g=2$ in our case) and $h$ is the Planck constant. It is important to stress that, the parameter $\epsilon_{c}$ serves to account for the finite size of galaxies. Note also that for $\epsilon_{c}\rightarrow+\infty$, we recover the Fermi-Dirac distribution. In the non-degenerate limit $\mu\rightarrow-\infty$, we recover on the other hand the classical King model King (1965), which reduces to the Boltzmann distribution in the limit $\epsilon_{c}\rightarrow+\infty$.

As we have mentioned, we are going to explore in this work the dynamical friction force effects on binary systems produced by DM profiles. However, it is important to note that, the dynamical friction force depends actually on the velocity distribution function whereby the introduction of the DM density profile, is somehow artificial; but in any case, it should be self-consistent for a given velocity distribution function according to the previous discussion.

III.3 The escape velocity

The escape velocity is defined in terms of the gravitational potential $\phi(r)$ of the background555Note that we are ignoring the gravitational potential produced by the binary system as well as other possibles contributions, as those produced by the baryonic component. as $v_{\rm esc}=\sqrt{-2\phi}$. The latter can be determined completely at any radius scale for a given density profile as follows

[TABLE]

The observed escape velocity of the Milky Way (considering the Galactic components, disk, bulge and halo) was found to be in the range $498$ km s*-1* $\lesssim v_{\rm esc}\lesssim 608$ km s*-1* at the solar position, at $90\%$ confidence interval and median likelihood of $544$ km s*-1* Smith et al. (2007). The RAVE survey has recently found the local escape speed to be $v_{\rm esc}=533^{+54}_{-41}$ km s*-1* Piffl et al. (2014). These values depend significantly on the mass exterior to the solar circle within a certain halo radius $r_{h}$. For example, the halo mass $M_{DM}(r_{h}=40$ kpc$)\sim 2\times 10^{11}M_{\odot}$ is consistent with the dynamics of the outer DM halo as was recently indicated in Gibbons et al. (2014). We note therefore that the Galactic escape velocity is either lower or closely equal to the orbital velocity of the binary pulsar for periods around $P_{b}\approx 0.1$ days. For large orbital periods $P_{b}\approx 100$ days, the orbital velocity is always well below the escape velocity. These two facts imply therefore that the contribution of the second integral (fast particles) to the dynamical friction force could be very small in most cases but not negligible in general. We will keep this term for a general study since, as we will see, it also leads to a change of sign in the orbital period time derivative (i.e. from decay to widening) for some values of the period as well as for the DM wind velocity.

III.4 The density profile

III.4.1 The NFW profile

We first recall the widely used phenomenological DM density profile arising within the $\Lambda$CDM cosmological paradigm, i.e. the NFW profile Navarro et al. (1996)

[TABLE]

where $\rho_{c}$ is the characteristic density and $r_{s}$ is the scale radius. This density profile exhibits a sharp cusp in the inner region $\rho\varpropto r^{-1}$ while in the halo part the density scales as $\rho\varpropto r^{-3}$.

It is worth to mention that there is an active debate in the literature on which is the best representation of the DM density profile that originates from the $\Lambda$CDM paradigm. For instance, some simulations has pointed out that the density profile of DM halos might be actually shallower than the one given by the NFW profile and found a cored structure represented more accurately by an Einasto profile (see Ref. Merritt et al. (2006) for details). It is out of the scope of this work to make an assessment on this issue and thus, for the sake of example, we adopt the NFW profile as the DM profile associated with the $\Lambda$CDM scenario. As we shall see, since the NSIS and the RAR profiles show also a shallower, cored inner halo666The similarity between the Einasto profile and the RAR profile in the inner halo region has been shown in Ref. Siutsou et al. (2015)., they are useful to analyze the differences that arise in the DMDF effect between cuspy and cored density profiles.

III.4.2 The NSIS profile

Another often adopted DM density profile which also yields the asymptotic flatness of the rotation curves is represented by the NSIS profile Shapiro et al. (1999):

[TABLE]

where $\rho_{0}$ is the central density and $r_{0}$ is the core radius.

III.4.3 The RAR profile

We will also examine the DMDF in the case of the RAR model Ruffini et al. (2015); Argüelles et al. (2016). This model describes the DM distribution along the entire galaxy in a continuous way, i.e. from the halo part to the Galactic center and without spoiling the baryonic component which dominates at intermediate scales. Likewise, the density $\rho$ and pressure $P$ for the Fermi-Dirac distribution function are defined respectively by

[TABLE]

Assuming a self-gravitating system of massive fermions (within the standard Fermi-Dirac phase-space distribution) in thermodynamic equilibrium, the DM density profile was computed in Ruffini et al. (2015). By imposing fixed boundary conditions at the halo and including the fulfillment of the rotation curves data, the parameters of the system have been constrained. This procedure was applied for different types of galaxies from dwarfs to big spirals exhibiting a universal compact core - diluted halo density profile. An extended version of the RAR model was recently presented Argüelles et al. (2016), by introducing a fermion energy cutoff $\epsilon_{c}$ in the fermion distribution. Importantly, this generalization in the statistics naturally arises by studying the stationary solution of a generalized Kramer statistics which includes the effects of escape of particles and violent relaxation Chavanis (2004). The new emerging density profile serves to account for the finite galaxy sizes due to the more realistic boundary conditions, while it opens the possibility to achieve a more compact solution for the quantum core working as a good alternative to the BH scenario in Sgr A* (see, Ref. Argüelles et al. (2016), for details). The narrow particle mass range provides several solutions to satisfy either the rotation curve data in the halo part or both sets of data, namely including additionally the orbits of the S-cluster stars such as the S2 star, necessary to establish the compactness of the DM central core. A comparison between the RAR model, NFW profile and NSIS for MW-like spiral galaxies is also shown in Fig. 2, describing the outstanding inner structure below parsec scale for the RAR profile.

It is important to clarify that the above DM density profiles are obtained without considering a DM-baryonic matter feedback nor DM self-annihilations. As shown in Ref. Bertone and Merritt (2005), these effects might produce changes in the DM density profile. We expect, however, the former to be important only locally in massive clusters and the latter stands on the largely model-dependent unknown DM nature. Thus, for the sake of generality, we shall not consider these effects in this work.

III.5 The velocity dispersion

According to observations of stars in outer part of halos and numerical simulation, the stellar velocity dispersion of the Milky-Way halo $\sigma_{r}$, shows an almost constant value around $120$ km s*-1* at scales of 20 kpc where DM is supposed to dominate and the circular velocity $V_{c}$ exhibits a flat behavior. Assuming that the galactic halo is stationary and spherically symmetric, it is possible to derive the DM radial velocity dispersion from the Jeans equation777Note that it is not the (observed) line of sight velocity dispersion of tracers.

[TABLE]

where $\beta=1-\sigma_{\theta}^{2}/\sigma_{r}^{2}$ is the velocity anisotropic parameter, that in the isotropic case, takes evidently the value $\beta=0$.

The circular velocity $v_{c}$ is defined by the local radial gradient of the potential while the radial velocity dispersion $\sigma(r)$ depends on the shape of the potential at exterior radii. For a non-rotating spherical system the relation between these quantities is given by

[TABLE]

where the first term in parenthesis is (minus) the logarithmic slope $\gamma$ of the density profile. For the singular isothermal sphere with a Maxwell Boltzmann distribution, the simple relation $v_{c}^{2}=2\sigma_{r}^{2}$ is satisfied for all radii. We note that, instead, for the NFW profile one obtains $1\leq\gamma\leq 3$ and hence, this simple relation between the circular velocity and the velocity dispersion is not fulfilled at all radii (except at the virial radius where $\gamma=2$, see e.g. Ref. Navarro et al. (2010)). Therefore, in order to find the right velocity dispersion profile for a given density profile, with associated gravitational potential, we solve hence the Jean equation for the isotropic case along the entire the Galaxy.

IV Orbital period evolution

In this section we study the DMDF effect as an intrinsic effect on the binary system motion. Hence, in order to analyze the perturbed Keplerian orbit of binary systems, we use the osculating formalism that permits to obtain the sequence of perturbed orbits Poisson and Will (2014). We follow particularly both the formulation and the derived analysis presented in Pani (2015) to compute the orbital period decay due to DMDF. We start by defining the relative acceleration between two bodies as

[TABLE]

with $\textbf{f}=a_{1}\eta\textbf{v}+a_{2}\textbf{v}_{w}$ for the case in which the perturbing force is taking to be the drag force measured on the center of mass. To zeroth order, the orbital velocity obeys a Keplerian motion, $v=\Omega_{0}r_{0}$, with $\Omega_{0}$ and $r_{0}$ being the angular velocity and orbital separation, respectively. We have also introduced the definitions: $\eta=\mu/M$, $\mu=m_{p}m_{c}/M$ and $M=m_{p}+m_{c}$. From here, the perturbed orbital elements can be then written as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the orbital parameters $a$, $e$, $\omega$, $i$ and $\Omega$ are the semiaxis major, the eccentricity, the longitude of pericenter, the inclination and longitude of the ascending node, respectively. In the right side of Eqs. (17)-(20), the source terms $S(t)$, $R(t)$ and $W(t)$ have been defined as a functions of the dynamical friction force as well as the wind velocity vector according to Pani (2015)

[TABLE]

[TABLE]

[TABLE]

The rate of change of the separation with time leads consequently to a change of the orbital period $P_{b}=2\pi/\Omega_{0}$ given by Landau and Lifshitz (1980)

[TABLE]

This relation along with Eq. (17) provide the time derivative of the orbital period888We note there is a typo in Eq. (18) of Ref. Pani (2015), namely when compared with Eq. (25) of our present work it shows an extra factor $v/2$ which leads the equation to be dimensionally incorrect.

[TABLE]

The resulting secular change in the orbital period is obtained by averaging over one period $P_{b}$, namely (see e.g. Ref. Pani (2015)):

[TABLE]

In the above formulation we have introduced the same definitions as in Pani (2015) for an easier comparison of the results: $\Gamma=v_{w}/v$, $\Delta_{\pm}=\Delta\pm 1$, $\Delta=\sqrt{1-4\eta}$. The coefficients $a_{i}$ can be written in terms of the integral velocity contribution function

[TABLE]

as

[TABLE]

with $A_{i}=4\pi\rho(r)\log\Lambda_{i}G^{2}M$. The definition of $b_{i}$ in the more general form given by Eq. (27) allows the use of any velocity distribution function, or equivalently any density profile through the integral term $I_{i}$ (term in parenthesis in Eq. (2)). This feature is contrary to the analyzed case in Pani (2015) where the Maxwell-Boltzmann distribution function was only considered there.

It is clear that the initial phase $\alpha$ can be set to any value without loss of generality, hence we set $\alpha=0$ for simplicity. In the next section V, we compute then the secular change of $\dot{P}_{b}$ for different density profiles with velocity dispersion profile determined by Eq. (14) and associated velocity distribution function described by Eqs. (7) and (8). The incorporation of the radial scale dependence of these quantities, leads to reduce the number of free parameters presented in early calculations Pani (2015), as already pointed out previously.

There may be other contributions to a secular change of the orbital period in addition to the DMDF and the gravitational-wave emission. A common effect in binaries with ordinary star components is the mass loss by star winds or accretion. A change of mass in the system would produce a change in the orbital period of the type $\dot{P}_{b}/P_{b}=-\dot{M}/M$, thus mass loss increases the orbital period (orbital widening) and mass accretion decreases it (orbital decay). In our case of binaries composed of compact stars the mass loss by winds is unlikely and accretion of matter from one component into the other could occur only via Roche lobe overflow for extremely short binary periods near the merging process. It remains the possibility of accretion of DM particles onto the binary components leading to a shrink of the orbit. The assessment of the importance of this effect, however, relies on the unknown cross-section between DM and baryonic matter inside the stars (see, e.g., Ref. Macedo et al. (2013)). Thus, for the sake of generality of our conclusions, we shall not include this effect in our estimates.

V Numerical results

We present now the dependence of $\dot{P}_{b}$ according to Eq. (25) on the free parameters: the orbital period, the DM wind velocity and the radial position of the binary measured from the Galactic center.

Once the density profile has been chosen, and the binary position has been fixed, the velocity distribution function, the velocity dispersion, as well as the escape velocity that constrains the maximum velocity in phase space, (upper limit in the second integral Eq. (2)) can be determined uniquely. Thus, for an observed binary at a known galactic position, the above quantities acquire values that can not be treated as uncorrelated and fully free parameters.

In the following analysis we adopt for the RAR model the solution for the Milky Way with a particle mass $m=345$ keV, which has the density profile with the most compact quantum core (see Fig. 2). We consider for the sake of example the following binary systems: NS-WD with masses $m_{p}=1.3M_{\odot}$ and $m_{c}=0.2M_{\odot}$, NS-NS with masses $m_{p}=m_{c}=1.3M_{\odot}$ and WD-WD with masses $m_{p}=0.5M_{\odot}$ and $m_{c}=0.25M_{\odot}$. According to our above discussion of the DM wind, and considering the observed orbital period range and binary positions, we perform our analysis varying the parameters in the following ranges: $10$ km s*-1* $\lesssim v_{w}\lesssim 1000$ km s*-1*, $0.1$ days $\lesssim P_{b}\lesssim 100$ days and a scale radius $0.1$ kpc $\lesssim r\lesssim 10$ kpc. It is important to note that we will also consider binary systems near the Galactic center (at parsec scales) since it is of interest to check the DMDF in regions along the Galaxy where DM is supposed to dominate.

V.1 DMDF in observed binaries

We first apply the approach to the Galactic binaries with measured intrinsic orbital periods and which are remarkably well explained by GW emission. In the last three columns of Table 1 we compare $\dot{P}^{\rm GW}_{b}$ with $\dot{P}^{\rm DF}_{b}$. In this calculation we use the NFW profile and the RAR model for illustrative purposes and the following free parameters: $\beta=0$ and $v_{w}=100$ km s*-1*. For other values of $\beta$, $\dot{P}_{b}^{\rm DF}$ does not change significantly, however a change of $v_{w}$ by one order of magnitude may be more important in the computation of $\dot{P}_{b}^{\rm DF}$ as we shall see below. At this point we should discuss whether the binaries of Table 1 are bound or unbound to the Galactic gravitational potential to determine a more precise value for the DM velocity wind101010It is important to clarify that the pulsar B2127+11C is located within the Galactic globular cluster M15 whereby it is subjected dominantly to the gravitational potential of its host globular cluster. As in the case of bulge globular clusters accelerating (possibly) pulsars through their stellar components Prager et al. (2016), DM can also contribute to the total acceleration by the studied effect in this paper. This latter claim is motivated by recent observational analysis that point out favorably the importance of the DM component in the dynamical of globular clusters Ibata et al. (2013); Sollima et al. (2016).. However, for the binaries of Table 1 which are characterized by short orbital periods, we checked that this is not relevant since for any DM wind in the range $10$ km s*-1* $\lesssim v_{w}\lesssim 1000$ km s*-1* the value of $\dot{P}_{b}^{\rm DF}$ is still very small compared with the $\dot{P}_{b}^{\rm GW}$ and with the measured intrinsic orbital period decay.

As we can see the DMDF effect is very small for all the above binaries because of the short orbital periods (compact orbits) that lead them to experience a small drag force.

We can thus first conclude that, for the binary systems listed in Table 1, the DMDF effect is indeed negligible and their secular evolution is fully dominated by GW emission.

V.2 DMDF as a function of the orbital period

A natural question that arises is whether DMDF effects can be comparable with the orbital period decay predicted by GW emission. To answer this question we explore the physical conditions (and hence the values for the model parameters) under which such equality may be attained. We thus consider the possibility to have binary systems with large periods, e.g. $P_{b}=100$ days, since DMDF is enhanced in systems with small binary compactness. We also consider regions along the Galaxy where the DM is supposed to dominate as those near the Galactic center.

We start our analysis by plotting the secular change of $P_{b}$ as a function of the orbital period for different density profiles in Fig. 3 with values for the free parameters $v_{w}=100$ km s*-1* and $r=0.1$ kpc in this analysis. We also show in the same plot the orbital decay due to GW emission $\dot{P}_{b}^{\rm GW}$, according to Eq. (1).

We can see that for a NS-WD system (top panel in Fig. 3), the orbital period decay starts to be dominated by the DMDF effect shortly after than $P_{b}=18$ days, for the NFW profile. i.e it is now larger that the one predicted by the GW emission. For the same system, The NSIS predicts a $\dot{P}_{b}^{DF}$ that matches $\dot{P}_{b}^{\rm GW}$ around $P_{b}=80$ days while for the RAR model, it occurs around $120$ days. For a NS-NS system (middle panel in Fig. 3), the NFW provides the match around $30$ days and around $90$ days and $130$ days for the NSIS and RAR model respectively. For a WD-WD system (bottom panel in Fig. 3), the NFW provides the match around $25$ days and around $70$ days and $150$ days for the NSIS and RAR model respectively. These results are also summarized in Table 2 for clarity. However, for such large periods, DMDF provides small orbital decays between $10^{-16}$ (for the NFW profile) and around $10^{-18}$ (for the other profiles) as can be seen in Fig. 3. It is evident from here that, the larger the orbital period, the larger the $\dot{P}_{b}$ reached. For instance for $P_{b}=1000$ days, $\dot{P}_{b}\sim 10^{-14}$ for the NFW profile and NS-NS binaries (middle panel in Fig. 3). These values are however very small, with respect, for instance, to the measured intrinsic orbital decays shown in Table 1 for some binary systems. However, possible measurements of the intrinsic period decays for binary systems with characteristic large periods is a challenge of unprecedented precision for astronomical observations. If such measurements might be successfully attained, it could also lead to discriminate between different DM density profiles due to the outstanding precision which is a characteristic property in such systems.

V.3 DMDF as a function of the DM wind

In order to analyze the effect of the wind velocity, we choose the radial position of the binary system fixed (measured from the Galactic center) at $r=1.5$ kpc and the orbital period $P_{b}=100$ days. The Fig. 4 shows that, for the aforementioned parameters and for the NFW and the NSIS profile, $\dot{P}^{\rm DF}_{b}$ lies in the range $10^{-20}$–$10^{-16}$. We can see from here that the smaller the DM wind velocity the larger the $\dot{P}^{\rm DF}_{b}$. However the latter statement does not apply for the RAR model which exhibits a constant value of $\dot{P}_{b}^{DF}\sim 5\times 10^{-18}$ for NS-WD and WD-WD and around $10^{-17}$ for NS-NS, for $v_{w}\gtrsim 200$ km s*-1*. this analysis leads to conclude that binaries into a DM background with small DM wind velocities (than the orbital velocity), experience a more effective drag force and hence a larger $\dot{P}_{b}$. We shall be then more interested in binary systems with small wind velocities, however we do not exclude at all binaries with (at least) one NS companion which may posse high kick velocities and then large wind velocities.

V.4 DMDF as a function of the binary position

We turn now to plot in Fig. 5 the value of $\dot{P}_{b}$ as a function of the radial position. We here adopt for the DM wind $v_{w}=200$ km s*-1* and for the binary period $P_{b}=100$ days (left panel) and $P_{b}=0.5$ days (right panel). We can see that, differences between the solution provided by $m=48$ keV and the one provided by $m=345$ keV for the RAR model is $\approx 3\times 10^{3}$. Interestingly, towards the Galactic center, for the two chosen cases of orbital periods ($P_{b}=0.5$, 100 days), the NFW profile and the RAR model (for the most compact solution $m=345$ keV) can reach a value of $\dot{P}_{b}$ that may be comparable with the one provided by Eq. (1) due to GW emissions. The prediction of $\dot{P}_{b}$ due to DMDF, for binaries with short orbital periods, can also be seen in Table. 1 for the NFW profile and the RAR model, respectively. For large periods however (left panel in Fig. 5), DMDF effect is highly enhanced for all the binary systems as can be seen in Fig. 3. In particular, the RAR model predicts large orbital period decay very near the Galactic center (around $10^{-3}$ pc) due to the high DM density at such distances (see also Fig. 2). The most promising situation arises then for binary positions near the Galactic center either for long or short orbital periods. We expect hence that observational measurements reach a technological improvement that permit us to measure such short orbital periods decays with outstanding precision in the future. In addition, it would be interesting to observe binary systems near the Galactic center to put constraints on the Galactic center environment, particularly on the DM density profile and importantly, to check the GR predictions in the strong field regime.

V.5 From orbital shrinking to widening

We turn now to analyze the model parameters under which a change of sign in the orbital period first time-derivative occurs. Namely, the conditions under which DMDF produces an orbital widening instead of an orbital shrinking or vice-versa. For given binary parameters and $\beta$, there are values of the wind velocity for which occurs a change of sign of $\dot{P}_{b}$. This is clearly seen in Fig. 4 for each density profile, for binaries with known values of the orbital period and distance and setting $\beta=\pi/2$. Fig. 6 shows, instead, how sensitive is this feature to the value of the $\beta$ parameter and to the Galactic DM distribution, i.e on the DM density profile. We describe now, for the sake of example and without loss of generality, the case of the NS-WD binary of Fig. 6. In this analysis we have adopted $P_{b}=100$ days and $r=1.5$ kpc as known quantities. The two changes of signs occur at: $\beta=68.75\degree$ and $114.6\degree$ for the RAR model with $v_{w}=70$ km s*-1*; $\beta=80.21\degree$ and $97.40\degree$ for the NFW profile with $v_{w}=200$ km s*-1*; and $\beta=74.49\degree$ and $103.13\degree$ for the NSIS profile with $v_{w}=300$ km s*-1* (see Fig. 6). These results are in general agreement with the ones found in Ref. Pani (2015) within the approximation of large $v_{w}$. The contribution of fast moving particles with respect to the binary-components, along with particular choices of $v_{w}$, $\beta$ and even $P_{b}$, might lead (although it is not a necessary condition) to multiple changes of sign of $\dot{P}_{b}$. This analysis supports the necessity of taking into account this contribution to check the conditions under which $\dot{P}_{b}$ may change sign. If one were interested in providing only negative values of $\dot{P}_{b}$, the particular choice $\beta=0$ (or more generally a value of it out the above ranges) would fulfill such requirement. We set then henceforth on the contrary $\beta=\pi/2$ in order to introduce a possible change of sign in $\dot{P}_{b}$ as a general case.

Let us turn back to Fig. 3. We note that the change of sign can occur for shorter or longer periods depending on the density profile and the DM wind velocity. For instance, for the NFW profile and $v_{w}=200$ km s*-1*, the change of sing occurs at $P_{b}\approx 2$ days contrary to the case $v_{w}=100$ km s*-1* where the change is around $100$ days for NS-WD as can be seen in the top panel of Fig. 3. The change of sign of $\dot{P}_{b}$ may occur at short period values for the RAR model, around 1 day for NS-NS and WD-WD (see middle and bottom panel in Fig. 3), while for the NFW profile, it may occurs around 20 days for NS-NS, around 12 days for WD-WD and around 100 days for NS-WD. The NSIS profile always provide negatives values in all the binary systems shown in Fig. 3. Before the first peak and after the second one, $\dot{P}_{b}$ is always negative while between the two peaks $\dot{P}_{b}$ is positive. We recall that, however, negatives values can be obtained for all the orbital period range in the case in which $\beta$ takes a value different from the aforementioned ranges, independently of the binary system and the other parameters as was inferred from Fig. 6.

It can be seen that negative values of $\dot{P}_{b}$ correspond to those binaries between the two peaks contrary to the curve given by the RAR model in the bottom panel in Fig 3. We also note that the position of those peaks does not change significantly when the orbital period varies, but rather the order of magnitude of $\dot{P}_{b}$. In some cases it can vary up to one order of magnitude. It is also important to note that this feature may change depending on the radial position and the density profile. As we can see from the same plot, the RAR model shows negative values of $\dot{P}_{b}$ below the first peak.

We can analyze the behavior as a function of the binary position in Fig. 5. For the case of NFW profile and left panel ($P_{b}=100$ days), negatives values of $\dot{P}_{b}$ can still be found after the peak as pointed out previously; therefore, positives values are located below $1.5$ kpc for the NFW profile, while for both the NSIS profile and the RAR model, $\dot{P}_{b}$ is always negative. In the right panel of the same figure (for $P_{b}=0.5$ days), all the DM density profiles provide negatives values of $\dot{P}_{b}$ except the RAR model, before the peak, for $m=48$ keV (also for $P_{b}=100$ days). It is important to stress that this analysis is valid for $\beta=\pi/2$ since, for other values of it, as $\beta=0$, $\dot{P}_{b}$ is always negative being independent of the DM density profile as can be inferred from Fig. 6.

VI Discussion and conclusions

It is by now well-known that the high-precision measurements of the orbital parameters of compact-star binaries (e.g. NS-NS, NS-WD and WD-WD) with short orbital periods ($P_{b}\lesssim 0.1$ days) have allowed a remarkable verification of the of the orbital decay predicted by general relativity due to GW emission (see Table 2 and references therein). However, the binary gravitational binding energy can be also affected by an usually neglected phenomenon, namely the DMDF (i.e. DM gravitational drag) induced by the DM on the binary owing to the interaction of the binary components with their DM gravitational wakes. We have qualified and quantified in this work this effect in the evolution of compact-star binaries and assessed the conditions under which it can become comparable to the one of the GW emission. We can draw the following conclusions from such an analysis:

A first interesting situation may occur for binaries with long orbital periods above 20 days: the orbital decay produced by DMDF becomes comparable to the one produced by the emission of GWs. Clearly, the precise orbital period at which the two effects are quantitatively equal depends on the DM density profile and on the binary parameters (see Fig. 3). 2. 2.

We have presented here, for the NFW, the NSIS and the RAR DM profiles, the orbital period for NS-NS, NS-WD and WD-WD binaries at which the DMDF effects, start to dominate over the produced by GW emission. These results are summarized in Table 2 (see also Fig. 3). 3. 3.

The NFW profile and the RAR model provide a more significant effect in the drag force than the one given by the NSIS profile, as can be seen in Figs. 3–5. It is important to note that the RAR and the NSIS profile predictions are similar above $P_{b}=100$ days, for all the binary systems analyzed in this work (with $v_{w}=100$ km s*-1* and located at $0.1$ kpc), while also for those values of $P_{b}$ the NFW profile predicts a much larger DMDF effect. 4. 4.

Another promising situation arises for binary systems located very near the Galactic center. In this case, the $\dot{P}_{b}$ due to DMDF is increased even for short orbital periods ($P_{b}=0.5$ days) as is shown in the right panel of Fig. 5. For long orbital periods the DMDF is notoriously strengthened, particularly for the NFW profile and the most compact solution for the RAR model ($m=345$ keV). This latter situation corresponds to the most ideal case for testing the DMDF (left panel of Fig. 5). 5. 5.

For the most ideal scenario of the DMDF effects in binary systems, kinematic effects, which are proportional to the orbital period, must be considered and respectively compared to the one studied in this work. 6. 6.

It is known that positive values of $\dot{P}_{b}$ can be caused for example by binary mass-loss or mass-exchange. However, we have seen that $\dot{P}_{b}$ might change sign from negative to positive due to DMDF. This is shown in Fig. 6 for different DM density profiles. Thus, this effect could be study in binary systems dominated by kinematic effects.

To summarize, The DMDF is very sensitive to the DM properties: density profile, velocity distribution function and velocity dispersion profile; whereby it would permit to put stringent constraints on the DM properties (and presumably on the nature) at the binary position and thus to discriminate between different DM models. Following this idea, the determination of the orbital secular changes of compact-star binaries with long/short orbital periods located in the outer halo/center of the Galaxy, might constrain the DM density distribution in these locations. It would be also interesting to study such an effect in binaries with measured orbital decays within globular clusters (as in the case of B2127+11C) in order to put constraints on the DM distribution in these systems. Therefore, the possible identification of this effect establishes a topic for future high-precision astrophysical data for the analysis of the secular evolution of compact-star binaries.

VII Acknowledgments

L.G.G. is supported by the Erasmus Mundus Joint Doctorate Program by Grant Number 2013–1471 from the EACEA of the European Commission. J.A.R. acknowledges the partial support of the project N 3101/GF4 IPC-11, and the target program F.0679 of the Ministry of Education and Science of the Republic of Kazakhstan. It is a pleasure to thank C. R. Arguëlles for interesting comments and discussions on the subject of this work.

Appendix A Distribution functions from the Eddington’s formula

For a given density profile, the gravitational potential can be obtained by solving the Poisson’s equation $\nabla\Psi=-4\pi\rho(r)$. Now, in order to solve the Eddington’s formula Eq. (6), we express the integral there in terms of $r$ instead of $\Psi$, and choose the appropriate limits of integrations by inverting numerically the equation $\Psi(r)=g(r)$, with $g(r)$ being a defined function of the radial position for a given density profile. In addition, the condition that the distribution function be positive for any positive energy, i.e, $f(\mathcal{E})\geq 0$ for $\mathcal{E}\geq 0$, should be guaranteed. This condition is fulfilled when $\Phi(r)$ goes to zero at infinity along with the appropriated value of the central potential $\Phi_{0}=\Phi(r=0)$.

We then set $\mathcal{E}=-E$ and do for convenience the simple change

[TABLE]

All quantities with bar are dimensionless by making use of the model parameters of the respective density profile. With all of this, we can perform numerically the integral in Eq. (6). In order to compare the distribution functions associated with the NFW and the NSIS profiles, with the Maxwell-Boltzmann one, we have to normalize $f$ to common units. For this we follow the usual normalization $\sqrt{8}M/(RV)^{3}$, where $M$ is the mass enclosed at a position $R$ where the circular velocity $V$ becomes flat, and $E$ is given in units of square velocity $V^{2}$; hence we introduce the dimensionless quantity $\bar{E}=E/V^{2}$. For the NFW profile (Eq. 10) such a radius is given by the virial radius, $r_{\rm v}=cr_{\rm s}$, where $c$ is the so-called concentration parameter $c$ and $r_{\rm s}$ is the scale radius. For this profile, we measure then $f$ in units of $\sqrt{8}M_{\rm v}/(R_{\rm v}V_{\rm v})^{3}$ and $\bar{E}=E/V_{\rm v}^{2}$. For the NSIS profile (Eq. 11), we adopt the core radius $r_{0}$ and thus all the quantities derived from it such that $f$ is given in units of $\sqrt{8}M_{0}/(r_{0}V_{0})^{3}$ and $\bar{E}=E/V_{0}^{2}$. Furthermore, we express the Maxwell-Boltzmann distribution as follows Binney and Tremaine (2008)

[TABLE]

with $\bar{E}=-\mathcal{E}/2\sigma^{2}$.

To check the consistency of our calculation we also apply the above method to the singular isothermal sphere (SIS)

[TABLE]

which must follow the the Maxwell-Boltzmann distribution function (see, e.g., Binney and Tremaine (2008)). We define the central density $\bar{\rho}_{0}=\sigma^{2}/2\pi G\bar{r}_{0}^{2}$ and compute its associated distribution function also from the Eddington’s formula. This solution is represented by the cyan-dotted line in Fig. 7 which can be seen overlaps with the Maxwell-Boltzmann distribution. For this profile, we measure $f$ in units of $\sqrt{8}\bar{M}_{0}/(\bar{r}_{0}\sqrt{2}\sigma)^{3}$. Thus, all the distribution functions and the dimensionless energy $\bar{E}$ are given in terms of theirs model parameters. Therefore, once we set the units of $f$ and $\bar{E}$, we can infer quantitatively the scale factor that lead to compare our results. However, it is important to mention that such a scale factor may be somewhat arbitrary when one does not consider a finite size for the halo which forces to introduce a cutoff at some radius scale. The relation between the relative energy $\mathcal{E}=-E$ and the particle velocity $v$ is determined by $E=\frac{1}{2}(v^{2}-v_{esc}^{2})$. For $v<v_{esc}$ particles are of course bounded. Finally, we present numerical results of the distribution functions associated to the NFW and the NSIS profiles and the comparison with the Maxwell-Boltzmann distribution function in Fig. 7. Our goal in this computation is to validate the approximation of taking the Maxwell-Boltzmann distribution to describe the velocity distribution for the aforementioned profiles. We can see that, the largest differences occur close to unbound energies, where precisely the contribution of particle velocities near the escape velocity do not contribute significantly to the dynamical friction force. These results then lead us to approximate, within our estimations, the velocity distributions function for the aforementioned profiles to follow the Maxwell-Boltzmann distribution. Such approximation permits to facilitate notoriously all the numerical computations regarding the orbital period decay. However, if we had at disposition observational timing pulsar data to test robustly our predictions, we would have to use the exact velocity distribution function for every density profile according to the Eddington’s formula.