Post-Newtonian Effects in Compact Binaries with a Dark Matter Spike: A Lagrangian Approach

Abstract

We apply the Lagrangian method to study the post-Newtonian evolution of a compact binary system with environmental effects, including a dark matter spike, and obtain the resulting gravitational wave emission. This formalism allows one to incorporate post-Newtonian effects up to any desired known order, as well as any other environmental effect around the binary, as long as their dissipation power or force formulae are known. In particular, in this work, we employ this method to study a black hole–black hole binary system of mass ratio ${10}^5$ by including post-Newtonian effects of order 1PN and 2.5PN, as well as the effect of relativistic dynamical friction. We obtain the modified orbits and the corresponding modified gravitational waveform. Finally, we contrast these modifications against the LISA sensitivity curve in frequency space and show that this observatory can detect the associated signals.

1. Introduction

The first direct detection of gravitational waves (GWs) by the LIGO/Virgo collaboration has opened up a new window into the universe [1]. Mergers of compact binary objects, such as black holes and neutron stars, provide unprecedented precision tests of general relativity and matter at its highest densities. There are also plans for space-based GW observatories such as The Laser Interferometer Space Observatory (LISA) [2], Taiji [3], and TianQuin [4]. These will be able to observe GWs at lower frequencies and, thus, to observe mergers of massive binary black holes and intermediate and extreme mass ratio inspirals (EMRIs/IMRIs). These systems are considered to be a rich source of signatures associated with various new and fundamental physics phenomena [5,6,7,8,9,10,11] and, in particular, dark matter (DM). The standard cosmological $\mathrm{\Lambda }$CDM model predicts the existence of such dark matter: a cold, collisionless massive particle that has so far eluded our direct detection efforts [12]. Utilizing this new window into the universe, GW probes of dark matter have been gaining traction in recent years [13].

When intermediate mass black holes (IMBHs) undergo adiabatic growth within DM halos, overdensities of DM, so-called DM spikes, can form [14,15]. We can find extensive work performed on modeling the spacetime of these structures, as well as their analytical forms given certain physical conditions, in Refs. [16,17,18,19]. The presence of these DM spikes can affect an inspiraling object as part of an IMRI. This was first explored in [20,21], where the DM spike interacts with the compact object through dynamical friction. This results in a faster inspiral compared to an inspiral in a vacuum and would be observable as a dephasing of the GW signal, possibly detectable by LISA [22,23].

Additional effects of these DM spikes have been explored in consecutive works, such as the effects of accretion of the DM spike [24,25], eccentric orbits inside these DM spikes [26,27,28,29], periastron precession [30], the halo feedback mechanism [23,31], relativistic corrections to dynamical friction and spike distribution [16,32], and DM spikes around primordial black holes [33]. These works all explore different effects that need to be combined in the end, as precise waveforms are needed to find them in the LISA data [34].

In this work, we present a general Lagrangian framework that can easily incorporate post-Newtonian (PN) corrections, dark matter dynamical friction, accretion, and any other orbital or environmental effect in a compact binary system with environment, as long as the mathematical formula of the aforementioned effects and corrections are known in the form of a dissipative power or a force. These effects are then formulated as generalized forces, and the Euler–Lagrange equations of motion yield the modified orbits, from which the waveform of the gravitational waves emitted by the binary system can be derived. Thus, this framework is different from other approaches that assume the quasi-adiabatic inspiral in their computations.

The structure of this paper is as follows. In Section 2, we present the system we are considering applying our Lagrangian formulation to and describe the model used for the DM spike profile. In Section 3, we present the main proposal of the paper, namely, the introduction of a Lagrangian formulation that can incorporate all the PN corrections up to any order, as well as other effects such as dynamical friction due to DM and mass accretion by the binary system. Section 4 is dedicated to studying the system detailed in Section 2, using the Lagrangian formalism of Section 4. There, we compute the orbits from which we derive the waveform of the GWs emitted by the compact binary system surrounded by a spike. We also compare the results with no-DM and no-PN cases and, by an analysis in the frequency space, show that these modifications to the GWs can be observed by LISA. Finally, in Section 5, we summarize our work and make some concluding remarks about the potential further applications of the Lagrangian framework we present in this paper.

2. Dark Matter Halo

2.1. Spike Profile

We consider a Schwarzschild black hole (BH) with mass $m_1$ that grows adiabatically and forms a surrounding DM spike ${\rho }_{\mathrm{DM}}\left(r\right)$ that concentrates DM from an initial Navarro–Frenk–White (NFW) profile [35]. A second smaller compact object of mass $m_2\ll m_1$ is on an approximately Keplerian orbit around the BH and experiences dissipative forces from gravitational wave emission and dynamical friction from the DM halo, which cause an inspiral due to the loss of orbital energy.

During the adiabatic growth of the central BH, a DM halo can contract and form a spike, resulting in large DM densities close to the BH horizon. The DM density profile describing such a spike was derived first in a semi-relativistic Newtonian manner [14] and later in a fully relativistic manner [15]. The fully relativistic model predicts that the DM density vanishes at $2R_s$ compared to $4R_s$ predicted by the semi-relativistic treatment, where $R_s$ is the Schwarzschild radius of the BH. Furthermore, the central densities of the relativistic DM spike can be significantly higher compared with the semi-relativistic case. Higher DM densities can have a significant impact on the rate of inspiral and, hence, the gravitational wave signal; therefore, we elect to use the fully relativistic model for our DM spike.

To model a relativistic DM density spike we follow the effective scaling function, Equation (7) in [32], given as

$${\rho }_{\mathrm{DM}}\left(r\right)=\overline{\rho }{10}^{\delta }{\left({\displaystyle \frac{{\rho }_0}{0.3\mathrm{GeV}/{\mathrm{cm}}^3}}\right)}^{\alpha }{\left({\displaystyle \frac{m_1}{{10}^6{\mathrm{M}}_{\odot }}}\right)}^{\beta }{\left({\displaystyle \frac{a}{20\mathrm{kpc}}}\right)}^{\gamma },$$

(1)

with

$$\overline{\rho }=A{\left(1-{\displaystyle \frac{4}{\tilde{x}}}\right)}^w{\left({\displaystyle \frac{4.17\times {10}^{11}}{\tilde{x}}}\right)}^q,$$

(2)

where $\alpha$, $\beta$, $\gamma$, $\delta$ are the relativistic NFW parameters generated by comparison with numerically generated curves; A, w, and q are fit parameters found by fitting to a reference curve with the scale parameters $({\rho }_0,M_{BH},a)=(0.3 \mathrm{GeV}/{\mathrm{cm}}^3,{10}^6{\mathrm{M}}_{\odot },20\mathrm{kpc})$; and $\tilde{x}=r/m_1$ (see Table 1 in Ref. [32])1.

For the inspiral process, the DM halo is assumed to be static. This assumption ignores the effects of dynamical friction on the DM halo itself, also referred to as halo feedback, which can be significant for binary systems with mass ratios less than ${10}^5$, as shown in Ref. [31]. Nevertheless, in this work, we first would like to check the consistency of applying the Lagrangian formulation to explore the effects of the post-Newtonian correction terms in the environment of a dark matter spike.

In the presence of a DM spike, the additional dissipative forces will speed up the inspiral, which is potentially observable in the GW signal. Therefore, it should, in principle, be possible to map the spike density. The application of match-filtering analyses such as the one in [33] are powerful tools for GW detections and parameter inference. It is, however, necessary to utilize suitable GW templates that, by taking into account a broader class of physical effects such as a relativistic DM spike or orbital evolution beyond the Newtonian regime, could improve the possibilities of potential detection. We consider the evolution of the binary system slightly before GW emission enters the lower end of the detector’s band. For our binary mass considerations, this coincides with regions of high density from the DM spike at a distance from the central BH of about a hundred times the innermost stable circular orbit, $r_{\mathrm{ISCO}}$, defined as

$$r_{\mathrm{ISCO}}=3R_{\mathrm{s}}.$$

(3)

2.2. System Parameters

In this work, we employ the static halo approximation, for which, following in the spirit of [31,32], we focus our study on a central mass of ${10}^6{M}_{\odot }$2 with a $10M_{\odot }$ companion. This produces a mass ratio of ${10}^5$, which is high enough such that a static halo model can be used as halo feedback, which becomes important at lower mass ratios, as well as ranging to the higher IMRI regime where most EMRIs are expected to have shed their dark matter halos [31]. To emphasize the sensitivity of the inspiral due to the DM density, we present three DM spikes in Section 4 with varying scale densities, ${\rho }_0$, i.e., ${\rho }_0=(0.1,0.3,0.5)\mathrm{GeV}/{\mathrm{cm}}^3$3. For simplicity, we match the scale radius to the reference value, i.e., $a=20\mathrm{kpc}$. Moreover, we assume a distance of $1\mathrm{Mpc}$ from the Earth to the binary and a total inspiral distance for the companion to be from $100r_{\mathrm{ISCO}}$ to $3r_{\mathrm{ISCO}}$, which is sufficient to show the effects of PN corrections, DM friction, while also being inside the range of LISA’s sensitivity band. Note that for these mass ranges (i.e., central black hole standing at ∼${10}^6{M}_{\odot }$), one may start to consider a self-force approach to model the waveform. However, we are considering the case where the evolution takes place at a relatively large distance of ∼$100r_{\mathrm{ISCO}}-3r_{\mathrm{ISCO}}$. With these conditions, the companion is situated far enough away from the strong field regimes and never exceeds $v/c\sim 0.3$, where PN correction terms are known to converge slowly [36]. An illustration of the dark matter density profile and the stages of the evolution for our system can be found on Figure 1.

3. Lagrangian Formulation and Equations of Motion

To incorporate post-Newtonian (PN) corrections and DM effects in our model to be able to compute the orbits with these corrections taken into account, we use the Lagrangian formulation and encapsulate PN and DM effects as generalized forces in this method. To see this more clearly, consider a system that is under both b Newtonian forces ${\tilde{\mathbf{F}}}^{\left(b\right)}$ and l non-Newtonian forces ${\mathbf{F}}^{\left(l\right)}$. The Lagrangian of the system can still be written as

$$L=T-V$$

(4)

where T is the kinetic energy of the system, and V is the potential energy associated to the Newtonian forces ${\tilde{\mathbf{F}}}^{\left(b\right)}$. However, using the d’Alembert principle, the Euler–Lagrange (E-L) equations of motion can be written as4

$${\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial {\dot{q}}^i}}\right)-{\displaystyle \frac{\partial L}{\partial q^i}}=\sum _l{Q}_i^{\left(l\right)}$$

(5)

where $Q_i^{\left(l\right)}$ represents generalized forces corresponding to non-Newtonian forces ${\mathbf{F}}^{\left(l\right)}$, i.e.,

$$Q_j^{\left(l\right)}={\mathbf{F}}^{\left(l\right)}·{\displaystyle \frac{\partial \mathbf{r}}{\partial q^j}}=F_i^{\left(l\right)}{\displaystyle \frac{\partial x^i}{\partial q^j}}.$$

(6)

Here, $\mathbf{r}$ is the Cartesian position vector of the object (or subsystem), and $q^j$ is its generalized coordinates, where in 3D, $j=1,2,3$. We can adopt $q^j$ to spherical coordinates such that $q^1=r,q^2=\theta ,q^3=\varphi$, and thus,

$$\mathbf{r}=\left(x^1,x^2,x^3\right)=\left(x,y,z\right)=\left(rcos\left(\varphi \right)sin\left(\theta \right),rsin\left(\varphi \right)sin\left(\theta \right),rcos\left(\theta \right)\right).$$

(7)

If the motion is restricted to a 2D plane, for example, in the case of a Keplerian orbit of a planet with $\theta ={\displaystyle \frac{\pi }{2}}$, then we can effectively work with just two coordinates $q^1=r,q^2=\varphi$, and thus,

$$\mathbf{r}=\left(x^1,x^2\right)=\left(x,y\right)=\left(rcos\left(\varphi \right),rsin\left(\varphi \right)\right).$$

(8)

For such a 2D system with l non-Newtonian forces, Equation (5) becomes

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial \dot{r}}}\right)-{\displaystyle \frac{\partial L}{\partial r}}=& \sum _l{Q}_r^{\left(l\right)},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial \dot{\varphi }}}\right)-{\displaystyle \frac{\partial L}{\partial \varphi }}=& \sum _l{Q}_{\varphi }^{\left(l\right)}.\hfill \end{array}$$

(10)

As mentioned above, in our model, we deal with both the PN corrections and the so-called dynamical friction effects. The latter stem from the interaction of the orbiting object (with lower mass) with its local DM environment and have nothing to do with friction in an electromagnetic sense.

In what follows, we add the conservative 1PN term and both dissipative 2.5PN and the dynamical friction terms, which are velocity-dependent, to the E-L equations of motion as generalized forces. Among the aforementioned PN terms, we first consider the leading terms that give corrections to the orbit and the energy dissipation through gravitational waves. The first orbit correction term is a conservative term related to the orbital precession, coming in at first order or 1PN [37,38]. The next leading term related to gravitational wave energy dissipation is a non-conservative one and is precisely equivalent to the so-called order 2.5PN correction term [37].

3.1. Newtonian Part of the Euler–Lagrange Equations

As mentioned above, in our two-body system that is surrounded by a DM spike, in the center of mass (CM) frame, we can write the Lagrangian as

$$L={\displaystyle \frac{1}{2}}\mu {\dot{\mathbf{r}}}^2-V\left(r\right)$$

(11)

where $\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2$ is the relative position of the objects in which indices $1,2$ refer to the central and the rotating objects, respectively, and ${\dot{\mathbf{r}}}^2=v={\dot{r}}^2+r^2{\dot{\varphi }}^2$. We denote the mass of the central more massive object by $m_1$ and the lighter orbiting object by $m_2$. Furthermore, the total mass of the two-body system is denoted by $m=m_1+m_2$ and its reduced mass by $\mu =m_1{m}_2/m$. The potential used in the Lagrangian is the Newtonian gravitational potential energy between the two bodies,

$$V\left(r\right)=-G{\displaystyle \frac{\mu m}{r}}.$$

(12)

Both masses can potentially be time-dependent if one also considers the accretion effect, by which they absorb DM from the DM spike around the central object.

3.2. Generalized Forces

Once we have the Lagrangian, we can write down the left-hand side of the E-L Equations (9) and (10). To write the right-hand side of these equations, we need to find the generalized forces. To carry that out, we assume that the corresponding usual forces themselves have two components: One component is (anti)parallel to the relative velocity of the objects $\mathbf{v}={\mathbf{v}}_1-{\mathbf{v}}_2$, and the second component is (anti)parallel to $\mathbf{r}={\mathbf{r}}_1-{\mathbf{r}}_2$, along the line connecting the two bodies. Such a force can be written in the form

$${\mathbf{F}}^{\left(l\right)}=F_r^{\left(l\right)}\widehat{\mathbf{r}}+F_v^{\left(l\right)}\widehat{\mathbf{v}}.$$

(13)

The corresponding generalized force is

$$Q_j^{\left(l\right)}={\mathbf{F}}^{\left(l\right)}·{\displaystyle \frac{\partial \mathbf{r}}{\partial q^j}}={\displaystyle \frac{F_r^{\left(l\right)}}{r}}\left(x{\displaystyle \frac{\partial x}{\partial q^j}}+y{\displaystyle \frac{\partial y}{\partial q^j}}\right)+{\displaystyle \frac{F_v^{\left(l\right)}}{v}}\left(\dot{x}{\displaystyle \frac{\partial \dot{x}}{\partial {\dot{q}}^j}}+\dot{y}{\displaystyle \frac{\partial \dot{y}}{\partial {\dot{q}}^j}}\right),$$

(14)

where we have used ${\displaystyle \frac{\partial \mathbf{r}}{\partial q^j}}={\displaystyle \frac{\partial \dot{\mathbf{r}}}{\partial {\dot{q}}^j}}$. Considering (8) and the fact that the index j in the above equation takes values $j=r,\varphi$, the radial and angular parts of the generalized force $Q_j^{\left(l\right)}$ above become (with $v=\sqrt{{\dot{r}}^2+r^2{\dot{\varphi }}^2}$)

$$Q_r^{\left(l\right)}=F_r^{\left(l\right)}+{\displaystyle \frac{F_v^{\left(l\right)}\dot{r}}{{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{\frac{1}{2}}}},$$

(15)

and

$$Q_{\varphi }^{\left(l\right)}={\displaystyle \frac{F_v^{\left(l\right)}{r}^2\dot{\varphi }}{{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{\frac{1}{2}}}}.$$

(16)

On the other hand, if we do not know the force itself but know the corresponding dissipative power $P^{\left(l\right)}$, we can write, using (13),

$$P^{\left(l\right)}={\displaystyle \frac{d{E}^{\left(l\right)}}{dt}}={\mathbf{F}}^{\left(l\right)}·\mathbf{v}=F_r^{\left(l\right)}\dot{r}+F_v^{\left(l\right)}{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{{\displaystyle \frac{1}{2}}}.$$

(17)

where we have used the polar coordinates expressions $\mathbf{v}=\dot{r}\widehat{\mathbf{r}}+r\dot{\varphi }\widehat{\mathit{\varphi }}$ and $\mathbf{r}=r\widehat{\mathbf{r}}$. In cases where the force is (anti)parallel to $\mathbf{v}$, i.e., when $F_r^{\left(l\right)}=0$, we obtain

$$F_v^{\left(l\right)}={\displaystyle \frac{P^{\left(l\right)}}{{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{\frac{1}{2}}}},$$

(18)

and (15) and (16) reduce to

$$\begin{array}{cc}\hfill Q_r^{\left(l\right)}=& {\displaystyle \frac{P^{\left(l\right)}\dot{r}}{{\dot{r}}^2+r^2{\dot{\varphi }}^2}},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill Q_{\varphi }^{\left(l\right)}=& {\displaystyle \frac{P^{\left(l\right)}{r}^2\dot{\varphi }}{{\dot{r}}^2+r^2{\dot{\varphi }}^2}}.\hfill \end{array}$$

(20)

3.3. Specific Form of Generalized Forces

As mentioned before, the generalized forces we consider in the model correspond to dynamical friction, $Q_j^{\left(\mathrm{DF}\right)}$, and two leading PN corrections, namely, the conservative 1PN correction to the orbit and the leading gravitational wave emission term, which is the 2.5PN correction $Q_j^{\left(\mathrm{GW}\right)}=Q_j^{\left(2.5\mathrm{PN}\right)}$.

3.3.1. Post-Newtonian Corrections

Following [37], we write the PN corrections to the relative acceleration of the two bodies perturbatively as

$${\displaystyle \frac{d\mathbf{v}}{dt}}={\displaystyle \frac{Gm}{r^2}}\left\{{\displaystyle \frac{1}{c^2}}{\mathbf{A}}_{1\mathrm{PN}}+{\displaystyle \frac{1}{c^4}}{\mathbf{A}}_{2\mathrm{PN}}+{\displaystyle \frac{1}{c^5}}{\mathbf{A}}_{2.5\mathrm{PN}}+{\displaystyle \frac{1}{c^6}}{\mathbf{A}}_{3\mathrm{PN}}+{\displaystyle \frac{1}{c^7}}{\mathbf{A}}_{3.5\mathrm{PN}}+\dots \right\},$$

(21)

with c being the speed of light in vacuum.

The 1PN term contribution above is of the form [37]

$${\mathbf{A}}_{1\mathrm{PN}}=\left\{(4+2\eta ){\displaystyle \frac{Gm}{r}}-(1+3\eta )v^2+{\displaystyle \frac{3}{2}}\eta {\dot{r}}^2\right\}\widehat{\mathbf{r}}+(4-2\eta )\dot{r}v\widehat{\mathbf{v}},$$

(22)

where $\eta =\mu /m$. Assuming the masses are time-independent, the corresponding 1PN force becomes

$$\begin{array}{cc}\hfill {\mathbf{F}}^{\left(1\mathrm{PN}\right)}=& \mu {\displaystyle \frac{d\mathbf{v}}{dt}}{|}_{1\mathrm{PN}}\hfill \\ \hfill =& {\displaystyle \frac{Gm\mu }{r^2{c}^2}}\left[(4+2\eta ){\displaystyle \frac{Gm}{r}}-(1+3\eta )v^2+{\displaystyle \frac{3}{2}}\eta {\dot{r}}^2\right]\widehat{\mathbf{r}}+{\displaystyle \frac{Gm\mu }{r^2{c}^2}}(4-2\eta )\dot{r}v\widehat{\mathbf{v}}.\hfill \end{array}$$

(23)

To find its corresponding generalized force, we replace (23) in (15) to obtain

$$Q_r^{\left(1\mathrm{PN}\right)}={\displaystyle \frac{Gm\mu }{r^2{c}^2}}\left[(4+2\eta ){\displaystyle \frac{Gm}{r}}-\left(3\eta +1\right)r^2{\dot{\varphi }}^2+{\dot{r}}^2\left(-{\displaystyle \frac{7}{2}}\eta +3\right)\right].$$

(24)

Likewise, from (23) and (16), we obtain

$$Q_{\varphi }^{\left(1\mathrm{PN}\right)}=(4-2\eta ){\displaystyle \frac{Gm\mu }{c^2}}\dot{r}\dot{\varphi }.$$

(25)

As stated before, the leading radiation term [38] is given by the 2.5PN correction, ${\mathbf{A}}_{2.5\mathrm{PN}}$. Again, assuming the masses are time-independent, the gravitational radiation term takes the form [37]

$${\mathbf{A}}_{2.5\mathrm{PN}}=-{\displaystyle \frac{8}{15}}\eta {\displaystyle \frac{Gm}{r}}\left\{\left(9v^2+17{\displaystyle \frac{Gm}{r}}\right)\dot{r}\widehat{\mathbf{r}}-\left(3v^3+9{\displaystyle \frac{Gm}{r}}v\right)\widehat{\mathbf{v}}\right\},$$

(26)

which results in

$$\begin{array}{cc}\hfill {\mathbf{F}}^{\left(2.5\mathrm{PN}\right)}=& \mu {\displaystyle \frac{d\mathbf{v}}{dt}}{|}_{2.5\mathrm{PN}}\hfill \\ \hfill =& -{\displaystyle \frac{8}{15}}\eta {\displaystyle \frac{G^{2}m\mu }{r^3{c}^5}}\left[\left(9v^2+17{\displaystyle \frac{Gm}{r}}\right)\dot{r}\widehat{\mathbf{r}}-\left(3v^3+9{\displaystyle \frac{Gm}{r}}v\right)\widehat{\mathbf{v}}\right].\hfill \end{array}$$

(27)

Using the above in (15) yields the radial part of the corresponding generalized force as

$$Q_r^{\left(2.5\mathrm{PN}\right)}=-{\displaystyle \frac{8}{15}}\eta {\displaystyle \frac{G^2{m}^2\mu }{r^3{c}^5}}\dot{r}\left[6\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)+8{\displaystyle \frac{Gm}{r}}\right].$$

(28)

In the same way, using (16), the angular part of this generalized force becomes

$$Q_{\varphi }^{\left(2.5\mathrm{PN}\right)}={\displaystyle \frac{8}{15}}\eta {\displaystyle \frac{G^2{m}^2\mu }{r{c}^5}}\dot{\varphi }\left(3\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)+9{\displaystyle \frac{Gm}{r}}\right).$$

(29)

3.3.2. Dissipation Due to Dynamical Friction

The gravitational interaction between $m_2$ and its local dark matter surrounding creates an effect known as dynamical friction (DF), similar in result to the electromagnetic friction, and slows down the rotating object, even though no electromagnetic or other forces in addition to gravity are assumed. This induces a dissipation of energy from the system, whose instantaneous power loss can be modeled as [21]

$$P^{\left(\mathrm{DF}\right)}={\displaystyle \frac{d{E}^{\left(\mathrm{DF}\right)}}{dt}}=-4\pi G^2{\displaystyle \frac{m_2^2{\rho }_{\mathrm{DM}}\left(r\right)}{v}}\xi \left(v\right)ln\left(\mathrm{\Lambda }\right).$$

(30)

Here, ${\rho }_{\mathrm{DM}}\left(r\right)$ is the dark matter density profile, v is the relative speed of the objects in polar coordinates, and $\xi \left(v\right)={\gamma }^2{(1+v^2/c^2)}^2$ is a relativistic correction to the DF in which $\gamma =\sqrt{1-v^2/c^2}$ is the familiar special relativistic factor. The term ${(1+v^2/c^2)}^2$ accounts for an increase in the deflection angle of the DM when considered as a collisional fluid due to the orbiting compact object, and ${\gamma }^2$ accounts for the relativistic momentum as seen by the compact object [39]. Furthermore, $\mathrm{\Lambda }$ is called the Coulomb logarithm and is defined as $\mathrm{\Lambda }=b_{\max }{v}_{\mathrm{typ}}^2/G\mu$, where $b_{\max }$ is the maximum impact parameter, and $v_{\mathrm{typ}}^2$ is the typical speed (squared) of the rotating object (or the relative speed squared in the CM frame). We take $\mathrm{\Lambda }=3$ [21].

Since this force is assumed to act only parallel to $\mathbf{v}$, we can use the above formula in (19) and (20) to obtain5

$$\begin{array}{cc}\hfill Q_r^{\left(\mathrm{DF}\right)}=& {\displaystyle \frac{P^{\left(\mathrm{DF}\right)}\dot{r}}{{\dot{r}}^2+r^2{\dot{\varphi }}^2}}=-4\pi G^2{m}_2^2\xi \left(v\right){\rho }_{\mathrm{DM}}\left(r\right){\displaystyle \frac{\dot{r}}{{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{\frac{3}{2}}}}ln\left(\mathrm{\Lambda }\right),\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill Q_{\varphi }^{\left(\mathrm{DF}\right)}=& {\displaystyle \frac{P^{\left(\mathrm{DF}\right)}{r}^2\dot{\varphi }}{{\dot{r}}^2+r^2{\dot{\varphi }}^2}}=-4\pi G^2{m}_2^2\xi \left(v\right){\rho }_{\mathrm{DM}}\left(r\right){\displaystyle \frac{r^2\dot{\varphi }}{{\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{\frac{3}{2}}}}ln\left(\mathrm{\Lambda }\right).\hfill \end{array}$$

(32)

3.4. Lagrangian Equations of Motion

We now have all the information to write down the Lagrangian equations of motion (9) and (10), which will take the form (with Newtonian parts on the left and non-Newtonian ones on the right)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial \dot{r}}}\right)-{\displaystyle \frac{\partial L}{\partial r}}=& Q_r^{\left(1\mathrm{PN}\right)}+Q_r^{\left(2.5\mathrm{PN}\right)}+Q_r^{\left(\mathrm{DF}\right)},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{\partial L}{\partial \dot{\varphi }}}\right)-{\displaystyle \frac{\partial L}{\partial \varphi }}=& Q_{\varphi }^{\left(1\mathrm{PN}\right)}+Q_{\varphi }^{\left(2.5\mathrm{PN}\right)}+Q_{\varphi }^{\left(\mathrm{DF}\right)}.\hfill \end{array}$$

(34)

Replacing (11), (12), (24), (25), (28), (29), (31), and (32) in the above yields

$$\begin{array}{cc}\hfill \ddot{r}-r{\dot{\varphi }}^2+{\displaystyle \frac{Gm}{r^2}}=& -{\displaystyle \frac{Gm}{c^2{r}^2}}\left[{\displaystyle \frac{16G\mu \dot{r}}{5c^{3}r}}+3\eta +1\right]\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)\hfill \\ \hfill & -{\displaystyle \frac{G\dot{r}}{c^2{r}^2}}\left({\displaystyle \frac{64G^2\mu m^2}{15c^3{r}^2}}-4m\dot{r}+{\displaystyle \frac{\mu \dot{r}}{2}}\right)+{\displaystyle \frac{2G^{2}m}{c^2{r}^3}}\left(2m+\mu \right)\hfill \\ \hfill & -{\displaystyle \frac{4\pi G^2{m}_2^2\xi {\rho }_{\mathrm{DM}}ln\left(\mathrm{\Lambda }\right)\dot{r}}{\mu {\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{3/2}}},\hfill \end{array}$$

(35)

and

$$\begin{array}{cc}\hfill r^2\ddot{\varphi }+2r\dot{r}\dot{\varphi }=& +{\displaystyle \frac{8G^2\mu m\dot{\varphi }}{5c^{5}r}}\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)+{\displaystyle \frac{2G\dot{r}\dot{\varphi }}{c^2}}\left(2m-\mu \right)\hfill \\ \hfill & +{\displaystyle \frac{24G^3\mu m^2\dot{\varphi }}{5c^5{r}^2}}-{\displaystyle \frac{4\pi G^2{m}_2^2\xi {\rho }_{\mathrm{DM}}ln\left(\mathrm{\Lambda }\right)r^2\dot{\varphi }}{\mu {\left({\dot{r}}^2+r^2{\dot{\varphi }}^2\right)}^{3/2}}}\hfill \end{array}$$

(36)

These coupled differential equations can be solved together numerically to yield the orbits.

4. Orbits and Gravitational Waves

4.1. Orbital Evolution

For our given parameters, to illustrate the effect of these PN corrections and DM friction, the semi-major axis $a\left(t\right)$ was computed with the numerical solutions to (35) and (36) as an evolution measure for all the effects introduced into the Lagrangian. These can be found in Figure 2.

As expected, DM friction accelerates the inspiral by several orders of magnitude. Note, however, that the addition of the 1PN term into the Lagrangian introduces a different orbital evolution, and in the presence of DM, it accelerates the orbital decay due to the effects of dynamical friction. This is particularly important during the later stages of the inspiral, where it differs from the regular GW radiative term. Note that the oscillations in the semi-major axis due to orbital precession from the 1PN term are still present at the later stages of the evolution; it is simply the scaling in plotting it alongside its no-1PN-corrections counterpart that makes the oscillations appear to vanish. Even though the 1PN term is a conservative term and would normally only introduce orbital precession, the presence of dark matter has affected the orbital evolution, which caused the orbit to lose energy faster than it would have with the normally used 2.5PN radiative term in the literature. One should expect that inclusion of higher order terms will also contribute to the inspiral’s evolution.

4.2. Gravitational Wave Analysis

Gravitational wave observatories measure the strain of the waves, which is the Fourier mode $h_{\sigma }$ of the perturbations $h_{ij}$, where $\sigma =+,\times$ are the polarization of the waves. Having found the orbits by solving the equations of motions (35) and (36), we replace the resulting $r\left(t\right)$ and $\varphi \left(t\right)$ in

$$\begin{array}{cc}\hfill x\left(t\right)=& r\left(t\right)cos\left[\varphi \left(t\right)\right],\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill y\left(t\right)=& r\left(t\right)sin\left[\varphi \left(t\right)\right],\hfill \end{array}$$

(38)

and use the Cartesian Dirac deltas to express the matter density of the system as

$$\rho (\mathbf{x},t)={\rho }_{\mathrm{DM}}\left(r\right)+\mu \delta (x-x\left(t\right))\delta (y-y\left(t\right))\delta \left(z\right).$$

(39)

Using these, we can find the quadrupole moment tensor

$$M^{ij}=\int x^i{x}^j\rho (t,\mathbf{x})d^{3}x,$$

(40)

where $i,j=1,2,3$ correspond to the $x,y,z$ coordinates, respectively. Finally, we can use $M^{ij}$ to compute the plus and cross polarizations of the gravitational wave strains using

$$\begin{array}{cc}\hfill h_{+}(t;\overline{\theta },\overline{\varphi })=& {\displaystyle \frac{1}{R}}{\displaystyle \frac{G}{c^4}}\left[{\ddot{M}}_{11}\left({cos}^2\left(\overline{\varphi }\right)-{sin}^2\left(\overline{\varphi }\right){cos}^2\left(\overline{\theta }\right)\right)\right.\hfill \\ \hfill & +{\ddot{M}}_{22}\left({sin}^2\left(\overline{\varphi }\right)-{cos}^2\left(\overline{\varphi }\right){cos}^2\left(\overline{\theta }\right)\right)\hfill \\ \hfill & -{\ddot{M}}_{33}{sin}^2\left(\overline{\theta }\right)-{\ddot{M}}_{12}sin\left(2\overline{\varphi }\right)\left(1+{cos}^2\left(\overline{\theta }\right)\right)\hfill \\ \hfill & \left.+{\ddot{M}}_{13}sin\left(\overline{\varphi }\right)sin\left(2\overline{\theta }\right)+{\ddot{M}}_{23}cos\left(\overline{\varphi }\right)sin\left(2\overline{\theta }\right)\right]\hfill \end{array}$$

(41)

and

$$\begin{array}{cc}\hfill h_{\times }(t;\overline{\theta },\overline{\varphi })=& {\displaystyle \frac{1}{R}}{\displaystyle \frac{G}{c^4}}\left[\left({\ddot{M}}_{11}-{\ddot{M}}_{22}\right)sin\left(2\overline{\varphi }\right)cos\left(\overline{\theta }\right)\right.\hfill \\ \hfill & +2{\ddot{M}}_{12}cos\left(2\overline{\varphi }\right)cos\left(\overline{\theta }\right)-2{\ddot{M}}_{13}cos\left(2\overline{\varphi }\right)sin\left(\overline{\theta }\right)\hfill \\ \hfill & \left.+2{\ddot{M}}_{23}sin\left(2\overline{\varphi }\right)sin\left(\overline{\theta }\right)\right],\hfill \end{array}$$

(42)

where R is the distance of the observatory to the CM of the binary system, and $\overline{\theta },\overline{\varphi }$ are related to the relative orientation of the frames of reference of the source and the observer [40]. Note that these are different from $\theta$ and $\varphi \left(t\right)$, which denote the angular position coordinates of the binary in the CM frame. A benefit to this approach is that by directly solving the Lagrangian equations of motion, one can freely compute the explicit waveforms as a function of time, where we can observe the rapid effect from DM friction through the dephasing of the waveforms. For our purposes, we have considered a source at $\{R,\overline{\theta },\overline{\varphi }\}=\{1\mathrm{Mpc},0,0\}$. Note that with this choice of angles, $h\left(t\right)=h_{+}\left(t\right)$. The result is plotted in Figure 3.

To fully encapsulate this GW dephasing effect, we analyze the number of cycles, ${\mathcal{N}}_{\mathrm{cycle}}$, that the system can be in-band as a function of the gravitational wave frequency, $f_{GW}$. Note that for this, we may use the relations that come directly from solving the Lagrangian equations of motion of the phase $\varphi \left(t\right)$, which is related to the GW phase with ${\varphi }_{GW}\left(t\right)=2\varphi \left(t\right)$, and the number of cycles would simply be

$$\begin{array}{c}\hfill {\mathcal{N}}_{\mathrm{cycle}}\left(t\right)\equiv {\displaystyle \frac{{\varphi }_{\mathrm{No}\mathrm{DM}}\left(t\right)-{\varphi }_{\mathrm{DM}}\left(t\right)}{\pi }}.\end{array}$$

(43)

Accompanied by the following relation,

$$\begin{array}{c}\hfill \dot{\varphi }\left(t\right)=2\pi f_{GW}\left(t\right),\end{array}$$

(44)

one may obtain ${\mathcal{N}}_{\mathrm{cycle}}\left(f_{GW}\right)$ with the help of the solutions to (35) and (36). We can see the result in Figure 4.

The number of cycles in the presence of a DM halo will be lower than in a vacuum, which is expected as the DM friction would translate a faster inspiral into fewer orbits, being in-band for GW emission. What is important for the detection of a DM spike from DM friction would be the particular shape of the waveform in frequency space, as interferometers set out to find signals from match-filtering techniques in the spectral domain. Hence, it is more illustrative to work in the frequency domain. In particular, we are interested in the detectability of these results with the LISA sensitivity curve, as well as the different signatures obtained from changing DM parameters. To adequately compare this to the detector’s sensitivity to estimate signal-to-noise ratios, we compute the characteristic strain $h_c\left(f_{GW}\right)$, given by

$$\begin{array}{c}\hfill h_c\left(f_{GW}\right)=2f_{GW}\left|\tilde{h}\left(f_{GW}\right)\right|\end{array}$$

(45)

with $\tilde{h}\left(f_{GW}\right)$ being the Fourier transform of the GW time signal. The numerical computation of the frequency domain in the case of the system we chose compared to the LISA sensitivity curve can be seen in Figure 5 (see also Appendix A). We can see there is a clear distinction between the cases with DM halo and the case without DM halo in lower frequencies, which will be observable by LISA.

To further accentuate these effects, looking at the ratio between the characteristic strain without any DM, or $h_{noDM}\left(f_{GW}\right)$, and the characteristic strain with dark matter, $h_{DM}\left(f_{GW}\right)$, we can see the large deviations mostly occur at the early stages of the inspiral, as shown in Figure 6. In particular, we are referring to the inspiral stages between $\sim 4·{10}^{-6}$ Hz to ${10}^{-4}$ Hz, or $100r_{ISCO}$ to $\sim 12r_{ISCO}$. This should be expected, as energy dissipation through GWs increases as the orbit shrinks. Eventually, GW emission dominates the energy loss from dynamical friction, and all profiles converge to the strain without dark matter. This, of course, means that the range of detectability of these DM environmental effects solely depend on the stage of the binary, where the signal would be maximal at lower frequencies, where clearly they will differ by at least an order of magnitude, as seen in Figure 6.

5. Discussion and Conclusions

Compact binary systems with an environment around them, particularly a DM spike, are a valuable experimental arena for the potential detection of effects associated with new physics, particularly the detection of the effects of DM on the emitted gravitational waves from such systems. Therefore, it is important to have a robust theoretical framework that can incorporate as many physical effects as possible, particularly general relativistic and DM effects. In this paper, we have employed a Lagrangian method that allows for the energetically consistent incorporation of conservative and dissipative effects, in particular, post-Newtonian corrections to the orbit up to any desired known order, as well as a relativistic dynamical friction term with DM. This allows us to expand the range of validity of the orbital evolution throughout the later stages of the inspiral. In particular, we applied this method to the case of a binary with a dark matter spike, where in addition to a relativistic form of dynamical friction, we have considered the leading PN corrections to the gravitational waveform, i.e., 1PN and 2.5PN terms, corresponding to the orbital precession and gravitational wave dissipation, respectively. After obtaining the analytical Lagrangian equations of motion, they were solved numerically to yield the orbits. These solutions were then used to compare the orbits in cases with and without DM halo and PN corrections. We have shown that the conservative 1PN term, while by itself will not lead to any orbital energy loss, in the presence of a DM spike, it will modify the evolution of the orbit, where in a wealth of literature, the dominant terms included are the GW dissipation term (equivalent to the 2.5PN term). Hence, any future realistic analysis should require corrections in order to improve the detectability of these effects. Furthermore, using the orbit equations, the modified GW waveform emitted from the binary in the presence of a DM spike was derived, which shows a dephasing effect compared to the no-DM case. Moreover, we have computed the frequency space strain of the GWs for three different DM densities and compared them against the LISA frequency sensitivity curve. Our calculations show that the modification effects resulting from the presence of a DM spike are observable by LISA.

In this method, the addition of further physical effects can be carried out in an energetically consistent manner through the Lagrangian and generalized forces. We have included a relativistic dynamical friction term and spike profile, as well as shown how to account for the addition of post-Newtonian corrections up to any known order and, in particular, corrections to both the orbital motion and GW emission up to the recently derived higher order terms [41]. Therefore, utilizing this framework could be quite beneficial in the search for DM effects, since following a potential inspiral beyond the Newtonian regime of the evolution will provide better GW templates for match filtering methods with the LISA detector or other observatories such as Cosmic Explorer. Moreover, this provides an alternative modeling tool to compare the current or future work in GWs sourced by DM [32,42]. One may also imagine modifying the Lagrangian to add other effects, such as the accretion of DM or baryonic matter by the binary, which can be incorporated into our framework quite easily by adding a time dependence on the masses. Lastly, this framework also provides the liberty to modify the environment density $\rho \left(r\right)$ to include dynamical and companion interaction effects $\rho \left(r\right)\to \rho (t,r,\varphi )$.