1 Introduction↩︎

One of the most intriguing open problems in the field of gravitational wave (GW) astronomy is the impact of so called “environmental effects” (EE) on GW signals. This term refers to the myriad of astrophysical and fundamental processes that can influence a GW source, alter its evolution and perturb its GW emission [1]–[28]. The study of EE is already considered a priority in the mission statement of the recently adopted Laser Interferometer Space Antenna (LISA) and the complementary proposed TianQin, since the majority of targeted GW sources, such as massive black hole (BH) binaries and extreme mass ratio in-spirals (EMRIs), are likely to merge deep within the accretion discs or nuclear clusters hosted in the centre of galaxies [29]–[35]. Additionally, recent work is strongly indicating that EE will leave detectable imprints in a significant fraction of stellar mass binary sources of GW for next generation ground based detectors [36].

The potential presence of EE in GW signals has profound implications on the science objectives of upcoming detectors, both in terms of fundamental physics and astrophysics. On the one hand, current waveform generation models are exclusively based on vacuum solutions of general relativity. If not accounted for, environmental influences can potentially confound results based on vacuum expectations, limiting the effectiveness of parameter estimation and introducing spurious biases, severely hinder the purported scientific goals of placing extremely stringent constraints on the strong field regime of gravity [7], [13], [37], [38]. On the other hand, environmental signatures in GWs can be used to probe the astrophysical surroundings of the sources they originate from. Prospects range from measuring properties of accretion discs [5], [17], [28], [39]–[41] or dark matter profiles [18], [42]–[45], to centre of mass accelerations [9], [46]–[52] and even tidal effects related to neutron star interiors [53]–[58]. Many of these measurements are beyond what is possible with other messengers, such as light or neutrinos, even in principle.

Most models of EE are based on the concept of dephasing, i.e. a phase shift between the vacuum signal and the signal including environmental perturbations. The standard procedure is to relate an environmental perturbation to a modification of the binary’s evolution in frequency, which in turn implies a change to the total accumulated phase of the GW signal. If the total dephasing is large enough, typically \(\sim \pi\) or \(\sim \pi/\)SNR [6], [21], [59], where SNR is the total signal-to-noise ratio of the signal, then the EE is said to be potentially detectable. Note that the magnitude of a detectable phase shift is also detector dependent. Most commonly, dephasing prescriptions reduce to power laws of the GW frequency \(f\), in a manner analogous to the various post-Newtonian (PN) corrections to the binary’s GW phase [60]–[62]. However, EE most typically induce dephasing with a much steeper frequency scaling, corresponding to the fact that they are more relevant for wider binary separations, where the intrinsic gravitational force of the binary is weaker. For this reason EE are said to appear at “negative PN orders", i.e. they accumulate with a steeper frequency scaling than even the leading order GW phase \(\propto f^{-5/3}\).

The detailed theory of building up and including PN corrections in GW templates -both in the time and the frequency domain- has been worked out in many seminal works in the 90s, for circular binaries [63]–[66]. Based on these works, the first investigations of EE have followed over the following decades, focusing in particular on neutron star (NS) tidal deformations [67]–[70] and the effect of gas drag on extreme mass ratio inspirals (EMRIs) [1], [5], [6]. Now, interest in the topic is rising due to the adoption of LISA and the advent of 3rd generation ground based detectors. More and more research on the astrophysical modelling side is demonstrating the relevance of EE in a plethora of GW sources, and more and more works are suggesting dephasing and other EE as an important observable to distinguish between binary formation channels and as a novel way to probe matter distributions in highly curved backgrounds [23], [36], [38], [49], [51], [71]–[74]. Furthermore, this trend is in concomitance with the increasing relevance of residual eccentricity in gravitational wave sources. Weak evidence for eccentricity is already present in a handful of LIGO-Virgo-KAGRA (LVK) signals [75], [76], and is considered a smoking gun for the dynamical and AGN formation channels for stellar mass binaries [77]–[89]. Similarly, EMRIs are most likely formed with high eccentricities [90], and recent research is showing evidence that even massive BH binaries may retain some eccentricity as they enter the mHz band [91]. All in all, the emerging picture is that a large quantity of signals in upcoming GW detectors will showcase complex features such as eccentricity, EE or most likely both. If the goal is to maximally benefit from such observations, it will be necessary to develop waveform templates capable of modelling these aspects to sufficient precision.

With this goal in mind, we present an overview of the fundamental concepts of how to include dephasing in GW waveforms, going from a description of EE in the time domain to a description of dephasing in the frequency domain, useful for GW data analysis. Our aim is to provide a comprehensive conceptual overview on the topic, as well as provide a useful reference for formulas and other prescriptions. We achieve this by revisiting and collecting the insight gained in many publications over the last two decades, as well as clear up some conceptual misunderstandings that are present in several more recent publications. Additionally, we lay some initial foundations for the analysis of dephasing due to EE in eccentric GW sources, which is a much richer problem with many open avenues of investigation. A depiction of the various aspects touched on by our work is presented in Fig. 1.

The paper is structured as follows. In section 2 we review how the Fourier domain waveform is obtained through the stationary phase approximation for circular binary inspirals, and interpret the connection between the time- and Fourier domain phases. In section 3 we provide the formulas and emphasize the important distinctions between different types of phase-shifts in the time- and Fourier domain. We also categorize the main types of environmental phase-shifts relevant for circular binaries, i.e. dephasing due to additional energy fluxes, radial potentials and Doppler shifts. In section 4 we go through some common astrophysical examples, in which several sources of dephasing act at the same time and compare their relative importance. In section 5 we extend the discussion of dephasing due to environmental effects to eccentric binaries, introduce novel ways in which energy and angular momentum fluxes can indirectly affect the GW phase, and illustrate the importance of these indirect effects through simple examples. In section 6 we close our overview with a brief discussion of signal-to-noise ratios related to the detectability of environmental dephasing. We finally summarize the main take-aways of our paper in section 7.

2 The stationary phase approximation - revisited↩︎

The analysis of GW signals can be performed in either the time domain -where the relevant quantity is the GW strain time series \(h(t)\)- or the frequency domain -where the relevant quantity is the Fourier transform of the strain \(\tilde{h}(f)\equiv{\rm FT}(h)\). While both representations are used in different applications, the numerical integration of GW signals is typically much more efficient in the frequency domain. In particular, the analysis of dephasing due to EE is almost exclusively performed in the frequency domain, by employing either simple estimates (see section 6), Fisher matrix methods [60], [92], or Bayesian inference methods [93]–[95]. The contrary is true when treating most perturbative effects on the binary inspiral. As an example, additional energy and angular momentum fluxes that modify the binary inspiral are typically given as functions of the binary parameters in time. Similarly, any other time-varying perturbations such as centre of mass accelerations, accretion, or other perturbative forces acting on the GW source are easy to model in the time domain. For these reasons, it is essential to be able to relate GW signals and dephasing prescriptions in the time domain to a representation of the perturbed signal in frequency space. The following ingredients are required: (i) the GW signal in the time domain (ii) a strategy to perform the Fourier transform analytically (iii) the relations between time and Fourier domain dephasing and (iv) a method to relate the dephasing to different physical scenarios. Let us first consider the first two points.

The detailed analysis of observed GW signals requires complex waveform models that contain all the relevant PN effects [96] and detector characteristics [97]. However, in order to assess the detectability of environmental effects it is often sufficient to incorporate dephasing prescriptions in the leading-order Newtonian GW waveform, neglecting here the problem of degeneracies. This is commonly done in phenomenological works, and is in fact more instructive as it allows to retain physical intuition about the GW signal and about the environmental dephasing itself. Note that this choice does not affect dephasing formulas in subsequent sections, where the relations \(t(f)\) and \(\phi(t_f)\) are not explicitly written out. Therefore, in this section we start by briefly summarising the formulas and the main takeaways of the seminal analysis of Newtonian frequency-domain GW waveforms, performed by Cutler & Flanagan [60]. The conclusions drawn are also valid for the general case. We use geometrical units with \(G=c=1\).

Let us also note that the dephasing prescriptions we provide here are valid when EEs introduce a small perturbation on the system, the evolution of which is driven mainly by GW radiation. It is also assumed here that the timescale of the frequency evolution, or "chirp" of the binary is short compared to the observation timescale. These conditions are not satisfied, for example, for wide binaries in triple systems, where the detectability of EEs requires a different treatment [19], [98], [99]. Additionally, other effects such as turbulent drag [100] or abrupt resonant excitations can instead introduce stochastic evolution or discontinuities, leading in the frequency domain to phase noise, spectral amplitude changes, or other distortions beyond a simple phase shift.

2.1 From time domain waveforms...↩︎

In the Newtonian description, binary BHs are considered as two point particles on a Keplerian orbit, which is secularly evolving due to quadrupolar GW radiation. The instantaneous orbital frequency is then \(\Omega=M^{1/2}/r^{3/2}\), where \(M = m_1 + m_2\) is the total mass. Then the evolution of the leading order quadrupolar GW frequency \(f\equiv f_\mathrm{GW} = \Omega/\pi\) reads: \[\frac{\mathrm{d}f}{\mathrm{d}t} = \frac{96}{5}\pi^{8/3} \mathcal{M}^{5/3}f^{11/3} \: , \label{eq:dfvac}\tag{1}\] where \(\mathcal{M} \equiv \mu^{3/5} M^{2/5}\) is the chirp mass and \(\mu = m_1 m_2/M\) is the reduced mass. Integrating this equation, we can obtain a closed solution for \(f(t)\). Then the observed GW strain becomes \[h(t) = \mathcal{A}(t)\cos\left(\int^t 2\pi f(t')\, \mathrm{d}t'\right) \: , \label{eq:ht}\tag{2}\] where \[\mathcal{A}(t)= \frac{(384/5)^{1/2}\pi^{2/3} Q\mu M}{D_Lr(t)}\] is the GW strain amplitude, which is a function of the component masses, the separation of the binary \(r(t)\), the luminosity distance \(D_L\) and the relative orientation of the binary and the detector, characterized by the scalar factor \(Q\) of \(\mathcal{O}(1)\). The prefactor in this equation is defined such that it is cancelled after a Fourier transformation.

Inverting \(f(t)\) we acquire the relation \(t(f)\) describing the time at which the GW signal reaches the frequency \(f\). This reads: \[t(f) = -5(8\pi f)^{-8/3} \mathcal{M}^{-5/3} + t_c \: , \label{eq:tf}\tag{3}\] where \(t_c\) is an integration constant having the meaning of a fictitious collision time, i.e. \(r(t_c) = 0\). In reality, the inspiral turns into a merger before the separation \(r\) can become formaly zero, producing a sharp cut-off in the GW signal. The time domain phase of the GW signal is given by the integral in the parenthesis of Eq. 2 . Evaluating the integral yields: \[\phi(t) = -2\left( \frac{t_c-t}{5\mathcal{M}} \right)^{5/8} + \phi_c \: ,\] where \(\phi_c\) is a second integration constant corresponding to the fictitious phase at the collision time. Expressing this as a function of the GW frequency we get yet another important relation: \[\label{eq:time95phase95freq} \phi(t(f)) = -2(8\pi\mathcal{M} f)^{-5/3} + \phi_c \:\tag{4}\] Eq. 4 denotes the time-domain phase, expressed as a function of frequency.

2.2 ...to Fourier domain waveforms↩︎

Now that we have the time domain waveform, we can relate it to its representation in the frequency domain via a Fourier transformation: \[\tilde{h}(f) \equiv \int\limits_{-\infty}^{\infty} e^{2\pi ift} h(t) \, \mathrm{d}t \: . \label{eq:FT}\tag{5}\] While the latter equation is formally correct, the integral cannot be evaluated analytically, even for the simplest circular Newtonian waveform. Therefore, some approximation scheme -the stationary phase approximation- is required to simplify this integral.

Circular binaries evolve slowly due to GW emission, with respect to the orbital timescale. Therefore, we can use that the modulation of the amplitude is slow with respect to the accumulation of the GW phase, i.e. \(\mathrm{d}(\ln \mathcal{A})/\mathrm{d}t \ll \mathrm{d}\phi/\mathrm{d}t\). Similar considerations apply to the GW frequency i.e. \(\mathrm{d}^2\phi/\mathrm{d}t^2 \ll (\mathrm{d}\phi/\mathrm{d}t)^2\). Under such conditions the integral can be evaluated using the stationary phase approximation, as follows: Given a time domain function \(B(t) = A(t) \exp[-i\phi(t)]\), we use the fact that the Fourier amplitude at a frequency \(f\) is mainly accumulated when the phase of the integrand in Eq. 5 becomes stationary, i.e. when \(\mathrm{d}\phi/\mathrm{d}t = 2\pi f\). This condition assures that the exponent in Eq. 5 is close to 1, such that the integral can accumulate rather than being dominated by the superposition of incoherent phases. Given the stationary phase condition, the phase in the integrand can be expanded as \[\label{eq:SPA95phase} 2\pi f t - \phi(t) \approx 2\pi f\, t_f - \phi_f - \frac{1}{2}\ddot{\phi}(t_f)(t-t_f)^2 \: ,\tag{6}\] where \(t_f \equiv t(f)\) and \(\phi_f \equiv \phi(t_f)\), and the linear term vanishes by construction. Then the integral can be evaluated analytically as a Gaussian integral, yielding: \[\begin{align} \tilde{B}(f) &\approx A(t_f) \left[ \frac{-i}{\dot{f}(t_f)} \right]^{1/2} \exp\left[ i(2\pi f t_f - \phi_f) \right] \nonumber \\ &= \frac{A(t_f)}{{\dot{f}(t_f)}^{1/2}} \exp\left[ i\left(2\pi f t_f - \phi_f - \frac{\pi}{4}\right) \right] \: . \end{align}\] The factor in front of the exponential can be understood intuitively. The Fourier domain amplitude of the waveform is proportional to the time domain amplitude, multiplied by the time the GW spends in that given frequency bin, which is proportional to \(\dot{f}(t_f)^{1/2}\). Thus the Fourier domain waveform in the Newtonian approximation becomes \[\tilde{h}(f) \approx \frac{Q}{D_L} \mathcal{M}^{5/6} f^{-7/6} \exp[i\psi(f)] \: ,\] where \(f \geq 0\) and \[\psi(f) = 2\pi f t(f) - \phi(t_f) -\frac{\pi}{4} \: . \label{eq:psif}\tag{7}\] We see that the Fourier phase \(\psi(f)\) is composed of two separate components. It contains information about the corresponding time domain GW phase \(\phi(t_f)\), but also about the time \(t(f)\) at which the binary reaches a given GW frequency. This will be a crucial insight when attempting to add perturbations of different types to the GW signal. Fig. 2 illustrates how we obtain the Fourier transform of the waveform in the stationary phase approximation. The waveform at a given time instance is approximated by a sinusoidal with the respective frequency, which has information about both \(t(f)\) and \(\phi(t_f)\).

2.3 General considerations↩︎

We now conclude our discussion of time and frequency domain waveforms with some additional important insights. While in the Newtonian example above we have given the exact functional forms of \(t(f)\) and \(\phi(t_f)\), Eq. 7 is true in general in the stationary phase approximation. Given a general prescription for \(\dot{f}\) we have \[\begin{align} \tag{8} t(f) &= \int^f \dot{f}^{-1} \, \mathrm{d}f' \: , \\ \tag{9} \phi(t_f) &= 2\pi\int^f f'\,\dot{f}^{-1} \, \mathrm{d}f' \: . \end{align}\] As \(\psi(f)\) contains two integration constants, let us take its derivative. Due to cancellations we have: \[\frac{\mathrm{d}\psi}{\mathrm{d}f} = 2\pi t(f) \: .\] Note that as \(\mathrm{d}\phi/\mathrm{d}t = 2\pi f(t)\), this makes \(t\) and \(f\) (or more precisely \(\omega=f/2\pi\)) conjugate variables and \(\psi(f)\) and \(\phi(t)\) related by Legendre transformations: \[\begin{align} \psi(f) &= \frac{\mathrm{d\phi}}{\mathrm{d}t}\bigg|_{t(f)} t(f) - \phi(t(f)) \: , \\ \phi(t) &= \frac{\mathrm{d\psi}}{\mathrm{d}f}\bigg|_{f(t)} f(t) - \psi(f(t)) \: . \end{align}\] This insight will be useful when interpreting the connection between dephasings in the time and Fourier domain.

Differentiating \(\psi(f)\) once again we obtain: \[\frac{\mathrm{d}^2\psi}{\mathrm{d}f^2} = \frac{2\pi}{\dot{f}} = \frac{8\pi^2}{\dot{E}_\mathrm{orb}} \frac{\mathrm{d}E_\mathrm{orb}}{\mathrm{d}\Omega}\: , \label{eq:dpsidf2}\tag{10}\] where in the second step we transformed the relation into a form used by multiple other studies [67], [101]. Crucially, Eq. 10 contains the full physical information about the Fourier phase evolution of the binary, while \(\psi(f)\) contains two additional degrees of freedom from the two integration constants.

Finally, we briefly note that the validity of the stationary phase approximation has been thoroughly studied in [102]. Corrections to the waveform appear in two forms: Firstly, the cubic term in the expansion used in Eq. 6 produces an additional phase term \(\delta\psi_{\rm SPA}^{\rm error}\) which is suppressed by 5PN orders: \[\begin{align} \delta\psi_{\rm SPA}^{\rm error} \propto \psi\times\left(\frac{v}{c} \right)^{10}. \end{align}\] Secondly, the finite observation time of real signals induces windowing effects that modify both the phase and the amplitude of the signal.

3 Dephasing in circular waveforms↩︎

We can now approach the topic of dephasing, i.e. of modelling perturbations on the GW signal related to some external influence acting on the binary inspiral. Firstly, note that deviations in the amplitude of the GW signal are extremely hard to produce, as they require a substantial instantaneous variation in the binary separation. On the contrary, even small instantaneous changes in the GW phase can accumulate over the entire inspiral and cause a significant difference in the overall waveform. There are many physical mechanisms that can induce phase shifts in the observed GW signals. In particular, perturbations to the orbital evolution induce a change in the binary’s chirp rate \(\dot{f}(f)\), which ultimately determines \(\psi(f)\). This can be realised as a result of additional energy fluxes that remove or inject energy from the orbit, or alternatively from the presence of external radial potentials, which change the binding energy of the orbit and the GW radiation at a given GW frequency. These two effects are usually called "dissipative" and "conservative" effects in GW modeling literature. Finally, observational effects, in particular Doppler shifts, can induce an apparent dephasing in the observed signal. In this section we discuss these different types of perturbations for circular inspirals, while also emphasizing and explaining the difference between time- and frequency-domain phase shifts.

The results of this section are summarised in Table 1, which lists the various relations between frequency chirp, time domain and Fourier domain dephasing for circular binary sources of GW.

3.1 Time-domain and Fourier domain dephasing↩︎

We start our discussion of dephasing with some general considerations that are independent of the cause of dephasing and the exact functional form of the perturbation on the binary. The only required simplification is for any perturbation to the binary inspiral to be small. Then, regardless of the specifics, we can express the binary’s total chirp rate as a sum of the vacuum chirp and a small perturbation \(\delta\dot{f}\): \[\dot{f}_\mathrm{tot}(f) = \dot{f}_\mathrm{vac}(f) + \delta \dot{f}(f) \: .\]

To obtain the dephasing in Fourier space we follow the same steps as in the vacuum case, but replace \(\dot{f}_\mathrm{vac}\) with \(\dot{f}_\mathrm{tot}\). Then \(t(f)\) and \(\phi(t_f)\) becomes: \[\begin{align} t(f) &= \int^f {\dot{f}_\mathrm{tot}}^{-1} \, \mathrm{d}f' = t_\mathrm{vac}(f) + \delta t(f) \: , \\ \phi(t_f) &= 2\pi \int^f \frac{f'}{\dot{f}_\mathrm{tot}} \, \mathrm{d}f' =\phi_\mathrm{vac}(t_f) + \delta\varphi(f) \: , \end{align}\] where \[\begin{align} \delta t(f) &\approx- \int^f\frac{\delta\dot{f}}{{\dot{f}_\mathrm{vac}}^2} \, \mathrm{d}f' \: , \\ \delta\varphi(f) &= \phi(t_f) - \phi_\mathrm{vac}(t_{f,\mathrm{vac}}) \approx - 2\pi \int^f\frac{f'\delta\dot{f}}{{\dot{f}_\mathrm{vac}}^2} \, \mathrm{d}f' \: . \end{align}\] Note that we have introduced the notation \(\delta \varphi(f)\) for the difference in the time-domain phases of the perturbed and vacuum signals at a given frequency to emphasize its distinction from the difference in the time-domain phases at a given time, \(\delta \phi(t)\). This makes the perturbation of the Fourier phase: \[\label{eq:dphasing95fourier} \delta\psi(f) = 2\pi f\delta t(f) - \delta\varphi(f) \: .\tag{11}\] Thus, the Fourier domain dephasing also contains two contributions, one from the binary dephasing at a fixed frequency, and one from the difference in time at which the binary reaches a given frequency. We note here that the additional contribution from the time shift \(\delta t(f)\) is not accounted for in multiple studies regarding dephasing in astrophysical scenarios [17], [21], [28], [103], [104]¹, or at least it is unclear methodologically whether it is [44], [72], [105], [106]. Nevertheless, the qualitative results of these works are most likely not affected, as the differences for power law dephasing prescriptions only account to pre-factors of order unity. We note that even though the time-domain phase shift contains all the necessary information about the EE, using the frequency-domain waveform usually simplifies the analysis, is more natural from a data-analysis perspective, and provides a better detectability criterion [107].

Differentiating Eq. 11 twice we obtain: \[\frac{\mathrm{d}^2(\delta\psi)}{\mathrm{d}f^2} = - 2\pi \frac{\delta \dot{f}}{\dot{f}_\mathrm{vac}^2} \: ,\] in analogy to Eq. 10 .

Figure 3 illustrates the various perturbed quantities, where we again highlight the difference between \(\delta\varphi(f)\) and \(\delta\phi(t)\). The perturbed orbit reaches the frequency \(f\) at a time \(t+\delta t\), meaning that: \[\begin{align} \delta\varphi(f) &= \phi_\mathrm{tot}(t+\delta t) - \phi_\mathrm{vac}(t) \nonumber \\ &\approx \phi_\mathrm{tot}(t) - \phi_\mathrm{vac}(t) + 2\pi f\delta t \: , \end{align}\] where we have assumed that \(\delta t\) is small. Thus: \[\delta\phi(t_f) \approx \delta\varphi(f) - 2\pi f \delta t(f) = - \delta\psi(f) \: . \label{eq:dphi95dpsi}\tag{12}\]

The above expression has important practical implications: It states that in order to obtain the magnitude of the Fourier domain phase shift of the GW signal it is sufficient to evaluate the time-domain phase difference at a specific time. This can expedite the evaluation of Fourier dephasing in binaries analysed in the time domain, as it allows to skip over a potentially costly Fourier transformation. In particular, such a relation may be useful in population synthesis codes, where large numbers of binaries are followed from formation to merger. Note also that since \(\phi(t)\) and \(\psi(f)\) are connected by a Legendre-transformation, Eq. 12 is analogous to saying that a Lagrangian perturbation \(\delta L\) generates a Hamiltonian perturbation \(-\delta H\).

Summarizing once again the different time and phase shifts for clarity:

• \(\delta t (f)\): time difference of the vacuum and perturbed signal reaching the reference frequency,

• \(\delta \varphi (f)\): time-domain phase difference between the vacuum and the perturbed signal at a given reference frequency \(f\),

• \(\delta \phi (t)\): time-domain phase difference between the vacuum and the perturbed signal at a given reference time \(t\),

• \(\delta \psi (f)\): frequency-domain phase difference between the vacuum and the perturbed signal at a given reference frequency \(f\).

3.2 Dephasing due to additional energy fluxes↩︎

We now start our discussion of individual mechanisms that induce dephasing in GW signals from inspiralling binaries, starting with the category of external perturbations. The simplest case is to directly change the binary inspiral rate by extracting or adding energy to it. In order to investigate these it is worth expressing \(\dot{f}\) in terms of the flux of orbital energy and the external perturbative energy flux \(\delta P\): \[\label{eq:enerflux} \dot{f}_\mathrm{tot} = \dot{E}_\mathrm{tot}\left(\frac{\mathrm{d}E_\mathrm{tot}}{\mathrm{d}f}\right)^{-1} = \dot{f}_\mathrm{vac} + \delta \dot{f}\tag{13}\] Therefore, a perturbation \(\delta P\) or \(\delta \dot{a}\) can be directly translated into \(\delta \dot{f}\): \[\begin{align} \label{eq:enerflux2} \dot{f}_\mathrm{tot} &= \left(\dot{E}_\mathrm{vac}+\delta P \right)\left(\frac{\mathrm{d}E_\mathrm{vac}}{\mathrm{d}f}\right)^{-1} \nonumber \\ &\approx\dot{f}_\mathrm{vac} \left[ 1 + \frac{\delta P}{\dot{E}_\mathrm{vac}} \right] \: , \end{align}\tag{14}\] thus: \[\delta \dot{f} \approx \dot{f}_\mathrm{vac}\frac{\delta P}{\dot{E}_\mathrm{vac}} \: . \label{eq:dfdP}\tag{15}\] Note however, that between Eqs. 13 and 14 –15 we have substituted the binary vacuum binding energy \(E_{\rm vac}\) for the binary total binding energy. This substitution is only valid in the case the external perturbation does not also induce a change in the binding energy, which is the topic of the following section.

3.3 Dephasing due to radial potentials↩︎

Beyond fluxes, the orbit of binaries can be perturbed by the presence of an external conservative potential \(\delta U\), as can be the case with three-body tidal interactions or the gravitational potential of a disc or dark matter spike. We show here that this will have multiple effects on the chirp of the signal. The clearest effect is the addition of an extra binding energy to the binary potential, which will have to be radiated away by GW emission. Additionally, the presence of an external conservative potential will modify the equilibrium point for circular binary orbits, i.e. binaries at a given radius will orbit with a different velocity in order to maintain a circular orbit. In turn, this induces a shift in the binary semi-major axis \(\delta a(f)\) that corresponds to a given GW frequency \(f\). The consequences are to modify the Keplerian binding energy, as well as implicitly change the GW radiation flux at a given GW frequency. In other words, one can understand this as the external perturbation causing a mismatch between the actual Keplerian binding energy or the GW radiation flux, and what an observer would expect by only considering the GW frequency \(f\). Thus, the relations \(E_\mathrm{vac}(a_\mathrm{vac}(f),f)\) and \(\dot{E}_\mathrm{vac}(a_\mathrm{vac}(f),f)\) no longer hold as the Keplerian relation \(a_\mathrm{vac}(f)\) does not hold. Here we consider the specific case of radial potentials \(\delta U(r)\). Then the total binding energy will read: \[E_\mathrm{tot}(f) \approx E_\mathrm{vac}(f) + \frac{\partial E_\mathrm{vac}}{\partial a}\Bigg|_f \delta a(f) + \delta U \: ,\] where \(E_\mathrm{vac}(f)\equiv E_\mathrm{vac}(a_\mathrm{vac}(f),f)\). The GW radiation flux will also become: \[\dot{E}_\mathrm{tot}(f) \approx \dot{E}_\mathrm{vac}(f) + \frac{\partial \dot{E}_\mathrm{vac}}{\partial a}\Bigg|_f \delta a(f) \: ,\] with \(\dot{E}_\mathrm{vac}(f)\equiv \dot{E}_\mathrm{vac}(a_\mathrm{vac}(f),f)\). This means that \[\begin{align} \frac{\delta \dot{f}}{\dot{f}_\mathrm{vac}} &\approx \dot{E}^{-1}_\mathrm{vac}\frac{\partial \dot{E}_\mathrm{vac}}{\partial a}\Bigg|_f \delta a(f) - \left( \frac{\mathrm{d}E_\mathrm{vac}}{\mathrm{d}f} \right)^{-1} \nonumber \\ &\times\left[ \frac{\mathrm{d(}\delta U)}{\mathrm{d}f} + \frac{\mathrm{d}}{\mathrm{d}f}\left(\frac{\partial E_\mathrm{vac}}{\partial a}\Bigg|_f \delta a(f)\right) \right] \: . \label{eq:delf95cons} \end{align}\tag{16}\] All of these different contributing terms must be considered when an external potential \(\delta U\) is present.

Let us try to simplify these relations by utilizing the Keplerian energy and 2.5PN GW radiation formulas. Since \(\dot{E}_\mathrm{vac}\propto a^4 f^6\), the first term becomes: \[\dot{E}^{-1}_0\frac{\partial \dot{E}_0}{\partial a}\Bigg|_f \delta a(f) = 4\frac{\delta a(f)}{a_0(f)} \: , \label{eq:delf95cons95rad}\tag{17}\] where we introduced the lower index \(0\) as shorthand for the vacuum definitions. From now on, for simplicity we also omit dependencies on \(f\) and \(a_0(f)\), where these are implied by the equations. Considering the force balance equations and defining \(\delta F \equiv-\mathrm{d(\delta U)/\mathrm{d}r}\) we get: \[\frac{\delta a}{a_0} \approx \frac{1}{6} \frac{\delta F}{E_0}a_0 \: .\] Let us now consider the second line in Eq. 16 , where we have: \[\frac{\mathrm{d(}\delta U)}{\mathrm{d}f} \approx -\delta F \frac{\mathrm{d}a_0}{\mathrm{d}f}\] For the partial derivative we have: \[\frac{\partial E_0(a,f)}{\partial a}\Bigg|_f = -4 \frac{E_0}{a_0} \: .\] Note the difference between this and the total derivative: \[\frac{\mathrm{d} E_0(a,f(a))}{\mathrm{d} a} = - \frac{E_0}{a_0} \: ,\] where the frequency is not fixed but is also changing according to Kepler’s law. Using these the second line in Eq. 16 , we continue the derivation: \[\begin{align} &\frac{\mathrm{d(}\delta U)}{\mathrm{d}f} + \frac{\mathrm{d}}{\mathrm{d}f}\left(\frac{\partial E_0}{\partial a}\Bigg|_f \delta a(f)\right) \nonumber \\ \approx& -\frac{5}{3} \delta F \frac{\mathrm{d}a_0}{\mathrm{d}f} - \frac{2}{3} a \frac{\mathrm{d(\delta F)}}{\mathrm{d}r} \frac{\mathrm{d} a_0}{\mathrm{d}f} \: . \end{align}\] Dividing this by \(\mathrm{d}E_\mathrm{vac}/\mathrm{d}f\) and further manipulating the equations we finally get: \[\begin{align} &- \left( \frac{\mathrm{d}E_0}{\mathrm{d}f} \right)^{-1} \left[ \frac{\mathrm{d(}\delta U)}{\mathrm{d}f} + \frac{\mathrm{d}}{\mathrm{d}f}\left(\frac{\partial E_0}{\partial a}\Bigg|_f \delta a(f)\right) \right] \nonumber \\ &\approx - \frac{10}{3} \frac{\delta U'}{V_g'}\Bigg|_{a_0(f)} + \frac{8}{3} \frac{\delta U''}{V_g''}\Bigg|_{a_0(f)} \: , \end{align}\] where \(V_g(r)\) is the gravitational potential and the prime denotes differentiation with respect to \(r\). After similar manipulations on Eq. 17 we get \[\dot{E}^{-1}_0\frac{\partial \dot{E}_0}{\partial a}\Bigg|_f \delta a \approx \frac{4}{3} \frac{\delta U'}{V_g'}\Bigg|_{a_0(f)} \: .\] Thus the total deviation in the binary’s frequency chirp becomes: \[\frac{\delta \dot{f}}{\dot{f}_\mathrm{vac}} \approx - 2 \frac{\delta U'}{V_g'}\Bigg|_{a_\mathrm{vac}(f)} + \frac{8}{3} \frac{\delta U''}{V_g''}\Bigg|_{a_\mathrm{vac}(f)} \: . \label{eq:Deph95U}\tag{18}\] It is worth further discussing the different terms in this equation. First of all, all the terms are proportional to only derivatives of the perturbing potential \(\delta U\). This is expected, as a constant shift in the potential should have no physical implications. The first term is the same as what we would get by only considering the effect of the extra binding energy from \(\delta U\). This is due to some accidental cancellations between the modification of the binding energy and the 2.5PN radiation power due to the fact that the orbit has a different semi-major axis at a given GW frequency. The second term, interestingly, is proportional to the second derivative of the perturbing potential and is due to the fact that \(\delta a\), which is proportional to the first derivative of the potential, is itself also changing as the orbit shrinks.

This equation has an additional interesting implication. We can make \(\delta \dot{f}/\dot{f}_\mathrm{vac}\) constant by introducing a potential with the same functional dependence on \(r\) as \(V_g\), i.e. \(\delta U\propto r^{-1}\). This, however, would be equivalent to increasing the total mass of the binary while keeping \(\mu\) constant. Note how the coefficient \(2/3\) we obtain from Eq. 18 in this case matches with the power-law dependence of \(\dot{f}_\mathrm{vac}\) on \(M\). On the other hand, the extra chirp can also be non-trivially cancelled up to first order by introducing a perturbing potential \(\delta U \propto r^{5/8}\). The dephasing due to e.g. a tidal potential linear in \(r\) would result in the same dephasing as what we would naively get from only taking into account the extra binding energy \(\delta U\).

3.4 Dephasing due to Doppler shift↩︎

The interpretation of the Fourier phase shift requires a slightly different treatment when considering a Doppler-shifted source. It is well-known that the time-domain strain, \(h(t)\), is invariant under the transformation [62]: \[(f,\mathcal{M},\mu,r,D_L,t) \rightarrow (f/\lambda,\mathcal{M}\lambda,\mu\lambda,r\lambda,D_L\lambda,t\lambda) \: .\] This means that a GW signal shifted by a constant Doppler shift \((1+z)\) would be indistinguishable from a non-Doppler-shifted signal with redshifted masses \(\mathcal{M}_\mathrm{obs} = (1+z)\mathcal{M}\) and \(\mu_\mathrm{obs} = (1+z)\mu\). However, time-dependent Doppler shifts (e.g. due to non-zero accelerations) cannot be modeled as a constant change in the binary parameter, but appear as frequency-dependent phase shifts, as we are about to show in this section. In order to get a more intuitive insight about this statement let us look at the Fourier phase of the observed signal of a binary source receding at a constant velocity \(v\). When the intrinsic GW frequency of the binary is \(f\), the observed frequency is Doppler-shifted by \[f_\mathrm{obs} = \frac{f}{1+v/c} \equiv \frac{f}{\lambda} \: .\] This part of the signal will also be observed at a different, Doppler-shifted time: \[t_\mathrm{obs}(f_\mathrm{obs}) = t(f) + \frac{v}{c}t(f) = \lambda t(f_\mathrm{obs}/\lambda) \: .\] Then, since for the time-domain phase \(\phi_\mathrm{obs}(t_\mathrm{obs}(f_\mathrm{obs})) = \phi(t(f))\), we have \[\begin{align} \psi_\mathrm{obs}(f_\mathrm{obs}) &= 2\pi f_\mathrm{obs} t_\mathrm{obs}(f_\mathrm{obs}) - \phi_\mathrm{obs}(t_\mathrm{obs}(f_\mathrm{obs})) \nonumber \\ &= 2\pi \lambda f_\mathrm{obs} t(f_\mathrm{obs}/\lambda) - \phi(t(f_\mathrm{obs}/\lambda)) \nonumber \\ &= \lambda^{-5/3} \psi^{\mathcal{M}}(f_\mathrm{obs}) = \psi^{\lambda\mathcal{M}}(f_\mathrm{obs}) \end{align}\]

Let us now relax the condition of constant velocity and investigate the phase shift for small perturbations. Then the binary signal with intrinsic frequency \(f\) is observed at a time: \[t_\mathrm{obs} = t(f) + \int^{t(f)} \frac{v(t')}{c}\mathrm{d}t' = t(f) +\Delta t(f)\] Then using that \(f=f_\mathrm{obs} + \Delta f\) we obtain the following for the Fourier phase: \[\begin{align} \psi&_\mathrm{obs}(f_\mathrm{obs}) = 2\pi f_\mathrm{obs} \left[t(f)+\Delta t(f)\right] - \phi(t(f)) \nonumber \\ \approx &2\pi f_\mathrm{obs}\left[ t(f_\mathrm{obs}) + \Delta f \frac{\mathrm{d}t(f)}{\mathrm{d}f}\Bigg|_{f_\mathrm{obs}} +\Delta t(f) \right] \nonumber\\ &- \phi(t(f_\mathrm{obs})) - \Delta f \frac{\mathrm{d}\phi(t(f))}{\mathrm{d}f}\Bigg|_{f_\mathrm{obs}} \: . \end{align}\] Using the definition of \(t(f)\) and \(\phi(f)\) in Eqs. 3 and 9 , we can see that the derivative terms cancel out. Then formally replacing \(f_\mathrm{obs}\) with \(f\) we get: \[\psi_\mathrm{obs}(f) \approx \psi(f) + 2\pi f\Delta t(f) \: ,\] where we used that \(\Delta t(f_\mathrm{obs}) \approx \Delta t(f)\) up to first order. Thus, we see that in contrast to our previous analysis of intrinsic dephasing, the Fourier domain dephasing of Doppler shifted GW sources is only due to the observational time difference between the vacuum and shifted binary, i.e. \[\delta \psi(f) = 2\pi f \Delta t(f) \: .\] This can also be shown by the equivalence relation in Eq. 12 . Since the time by which the observed signal reaches the phase \(\phi(f)\) is delayed by \(\Delta t(f)\), the phase of the Doppler-shifted signal lags behind by: \[\delta \phi(t_f) \approx -2\pi f_\mathrm{obs} \Delta t(f) \approx -2\pi f \Delta t(f) \: ,\] which agrees with our previous result.

Also note that since a constant velocity simply appears as a change in the chirp mass (as shown in the beginning of the section), the actual observable dephasing would be due to a deviation from this constant velocity, i.e. an acceleration of the center of mass. Denoting our reference time for the constant recession as \(t_0\), we have: \[\Delta t(f) = \int\limits_{t_0}^{t(f)} \frac{v(t')-v(t_0)}{c} \mathrm{d}t'\]

This formula is also applicable to redshifts not originating from proper motion, i.e. gravitational or cosmological redshift. In this case the time difference translates to a difference in proper and observer time and the formula becomes: \[\Delta t(f) = \int\limits_{t_0}^{t(f)} \left[ z(t')-z(t_0) \right] \mathrm{d}t'.\]

4 Illustrative Examples↩︎

In realistic astrophysical environments phase shifts typically originate from multiple sources. For example, perturbations from a third body will induce a center-of-mass acceleration, create an additional conservative perturbing potential, and induce secular changes in the orbital parameters. Drag forces from the surrounding gas will act as a source of extra energy flux, induce an additional conservative potential, while asymmetric forces on the two components will also accelerate the center of mass of the binary [108].

We now go through some simple analytic examples in astrophysics, in which several sources of dephasing act at the same time. We split the discussion into two separate astrophysical sources of GW, i.e. extreme mas ratio inspirals (EMRIs) in the first part and stellar mass binary sources in the second. For the former, we analyse the relative importance of drag vs. conservative radial potentials in typical matter backgrounds. For the latter we compare the importance of Doppler shifts due to accelerations vs. gravitational redshifts in the presence of a third massive body.

4.1 Drag vs potentials in EMRI inspirals↩︎

EE in EMRI sources are the subject of many studies [5], [7], [13], [17], [18], [28], [37]–[40], [42]–[45]. The two most commonly studied effects are 1) the presence of gas torques for EMRIs embedded in accretion discs, and 2) the effect of dynamical friction in dark matter spikes. Here, we compare and discuss the dephasing caused by the two aforementioned astrophysical effects with the dephasing caused by the additional potential \(\delta U\) implied by the presence of a gaseous disc or a dark matter distribution (also called spike) around the central massive BH.

Accretion disc potentials and dark matter distributions around massive black holes are a field of active research. Here we employ two simplified models, that capture the basic physics at play. As a proxy for a realistic accretion disc, we use the potential and surface density profile of a rigidly rotating Mestel disk [109]–[111]. The surface density \(\Sigma\) of a Mestel disk is given by the simple relation: \[\begin{align} \Sigma = \frac{M_{\rm d}}{2 \pi R_{\rm d}}\frac{1}{r} \end{align}\] where \(M_{\rm d}\) is the disc’s total mass, \(R_{\rm d}\) its truncation radius and \(r\) is a radial coordinate. The potential in the midplane is given by: \[\begin{align} \delta U_{\rm d}(R) = -\frac{GM_{\rm d}}{R_{\rm d}} \ln\left(r/R_{\rm d} \right). \end{align}\] For sufficiently large radii, the potential experienced by the EMRI is therefore well approximated by the Newtonian result: \[\begin{align} U_{\rm tot} \approx -\frac{G M_{\bullet}}{r} -\frac{GM_{\rm d}}{R_{\rm d}} \ln\left(r/R_{\rm d} \right), \end{align}\] where \(M_{\bullet}\) is the central BH mass. As shown before, the presence of this additional potential would induce a perturbation in the frequency chirp. We evaluate Eq. 18 to find:

\[\begin{align} \delta \dot{f}_{\rm pot} \approx q\frac{64 \pi^2}{5 c^5}\frac{G^3 M_{\bullet} M_{\rm d}}{R_{\rm d}}f^{3}, \end{align}\] where \(q <<1\) is the EMRI mass ratio. This additional chirp contribution increases rapidly as a function of frequency. This corresponds to the fact that, as the EMRI chirps faster and faster radially due to GW radiation, it is experiencing more variation in the perturbative potential as a function of time.

We now compare the dephasing due to the perturbative potential to the effect of gas torques. Small orbiters embedded in thin accretion discs excite density perturbations, which can induce a back reaction on the mass. This effect is known as migration torque, or type I torque, and is commonly invoked as an important EE. The magnitude of type I migration torques is [112]: \[\begin{align} \Gamma_{\rm I} = -q^2\Sigma r^4 \Omega_{\rm K}^2 \mathcal{M}_{\rm a}^2 \end{align}\] where here \(\Omega_{\rm K}\) is the Keplerian angular frequency and \(\mathcal{M}_{\rm a}\) is the disc’s Mach number. Here we assume a constant Mach number, which corresponds to isothermality for a rigidly rotating Mestel disc. Type I torques induce an additional energy flux onto a circular orbit: \[\begin{align} \delta P_{\rm I}= \Gamma_{\rm I} \Omega_{\rm K}. \end{align}\] As shown in Eq. 15 , this implies a corresponding change in the frequency chirp: \[\begin{align} \delta \dot{f}_{\rm I}= 3q \frac{\mathcal{M}_{\rm a}^2 G^{2/3}\Sigma}{\pi ^{1/3} M_\bullet^{1/3}} f^{2/3}. \end{align}\] Evaluated for our simple disc model, we find: \[\begin{align} \delta \dot{f}_{\rm I}= q \frac{3\mathcal{M}_{\rm a}^2 G^{1/3} M_{\rm d}}{2\pi ^{2/3} M_\bullet^{2/3} R_{\rm d}} f^{4/3}. \end{align}\]

The ratio \(\mathcal{R}_{\rm d}\) between the induced frequency chirps from a Mestel disc reveals some insight: \[\begin{align} \mathcal{R}_{\rm d}= \frac{\delta \dot{f}_{\rm pot}}{\delta \dot{f}_{\rm I}} = \frac{128\pi^{8/3}}{15 \mathcal{M}_{\rm a}^2} \left(\frac{G M_{\bullet}}{c^3} f\right)^{5/3}, \end{align}\] which rescaled in terms of the ISCO frequency for a non-rotating BH is: \[\begin{align} \mathcal{R}_{\rm d}= \frac{1}{\mathcal{M_{\rm a}}^2} \frac{16 \pi\sqrt{2} }{135\sqrt{3}}\left(\frac{f}{f_{\rm ISCO}} \right)^{5/3}. \end{align}\] We see that even for highly relativistic EMRIs, the dephasing caused by the disc’s potential is suppressed with respect to the gas torque effect by a factor \(\mathcal{M}_{\rm a}^{-2}\). Thus, for any but the hottest of discs it is entirely neglectable, and we recover the conclusions discussed in [6], [7].

We now turn our attention to the case of dark matter spikes. Here we model the density of the spike as a singular isothermal sphere: \[\begin{align} \rho = \rho_0 \left(\frac{r}{r_{\rm ISCO}}\right)^{-2} \end{align}\] where \(\rho_0\) is typically \(10^{18}\) GeV/cm\(^3\) to \(10^{20}\) GeV/cm\(^3\), scaled to the radius of the ISCO [44]. The enclosed mass of the spike is: \[\begin{align} M_{\rm encl} = r\times 4 \pi \rho_0 r_{\rm ISCO}^2 \end{align}\] which implies a gravitational potential: \[\begin{align} \delta U_{\rm DM} = -4 \pi G\rho_0 r_{\rm ISCO}^2 \ln(r/r_{\rm ISCO}). \end{align}\] As we can see, this has an identical functional form as the Mestel disc we previously discussed. In dark matter spikes, a small orbiter experiences a drag force due to dynamical friction. In the supersonic limit, dynamical friction is well approximated by the following formula [113]: \[\begin{align} \lvert \mathbf{F}_{\rm dyn} \rvert = q^2\frac{4 \pi G^2M_{\bullet}^2}{v^2}\rho\ln(\Lambda), \end{align}\] where \(v\) is the orbital velocity and \(\ln(\Lambda)\) is the Coulomb logarithm. This induces a modification to the EMRI chirp, in complete analogy to the disc embedded case. The energy flux is \(-\lvert \mathbf{F}_{\rm dyn} \rvert v\) and the frequency chirp is: \[\begin{align} \delta\dot{f}_{\rm DF}= 12\, q\, G \rho(f) \ln(\Lambda). \end{align}\] We compare this with the potential induced chirp, which we compute by evaluating Eq. 18 :

\[\begin{align} \delta\dot{f}_{\rm DM} = - 2304\pi^3 \frac{G^2 M_\bullet\rho_0 r_{\rm ISCO}^2}{5 c^5}qf^3. \end{align}\] Once again, the ratio \(\mathcal{R}_{\rm DM}\) of the frequency chirps induced by the presence of a dark matter spike is informative: \[\begin{align} \mathcal{R}_{\rm DM} = \frac{\delta \dot{f}_{\rm DM}}{\delta \dot{f}_{\rm DF}} &= 192 \pi^{5/3}\left( \frac{G M_{\bullet}}{c^3} f\right)^{5/3} \nonumber \\ &= \frac{8}{15}\sqrt{\frac{2}{3}} \left( \frac{f}{f_{\rm ISCO}}\right)^{5/3}. \end{align}\] where we neglect the Coulomb logarithm. We see here that for highly relativistic EMRIs, both contributions to the dephasing can be of comparable order. This implies that the two contributions to dephasing in such systems must be inferred jointly, in order to properly recover the value of the dark matter density.

Note how the conclusions for Mestel discs and isothermal sphere dark matter mostly differ by a factor of the Mach number squared, which is purely a measure of the geometrical flattening of the matter distribution. This hints at a broader connection between dynamical friction – caused by the gravitational pull of wake-like overdensities – and type I torques, caused by the gravitational pull of interference patterns between wake-like overdensities.

4.2 Doppler and gravitational redshifts↩︎

Here we compare two types of redshifts that can modify the observed frequency of GW sources. As discussed previously, dephasing due a Doppler shift of the GW signal can be observed if the inspiralling binary is accelerated along the line of sight. We consider here the acceleration caused by a third massive object with mass \(m_3\) and at a distance \(r\). The change in redshift caused by the line-of-sight acceleration \(a_{\rm ls}\) is: \[\dot{z}_{\rm D} = \frac{a_{\rm ls}}{c} = \frac{Gm_3 }{r^2c} \times P \: ,\] where \(P\) is the projection factor onto the line of sight and we assume that the acceleration is constant. Thus the phase shift is: \[\delta \psi_\mathrm{D} = 2\pi f \frac{Gm_3}{r^2c} \frac{t^2}{2} P \: .\] Using Eq. 3 and \(c=G=1\) we get at leading order (cf. Eq. (9) in [46] or Eq. (4) in [47]) \[\delta \psi_\mathrm{D} = \frac{25 M}{65536 \eta^2} a_\mathrm{ls} (\pi M f)^{-13/3} \: ,\] where \(\eta=\mu/M\) is the symmetric mass ratio.

Thus, a Doppler induced dephasing arises as a consequence of the acceleration induced by the third body. However, the presence of a third body has the additional implication of changing the gravitational potential in the vicinity of the GW source, inducing a gravitational redshift [9], [20], [114]. Taking into account the time dilation effect due to the proper motion of the center of mass as well, we get the following expression for weak gravitational fields [114]: \[\begin{align} 1+z_\mathrm{G} = \left( -g_{\mu\nu} V^\mu V^\nu \right)^{-1/2} \approx 1 + \frac{Gm_3}{c^2 r} + \frac{V^2}{2c^2} \: , \end{align}\] where \(g_{\mu\nu}\) is the Schwarzschild metric [115], \(V^\mu = (1,\boldsymbol{V})\) with \(\boldsymbol{V}\) being the center-of-mass velocity of the binary, and \(r\) is the distance from the third object with mass \(m_3\). Then for a binary moving radially: \[\begin{align} \dot{z}_\mathrm{G} \approx - \frac{Gm_3}{c^2r^2} V_r + \frac{a_r V_r}{c^2} = -2\frac{Gm_3}{c^2r^2} V_r \: , \end{align}\] where \(V_r\) is the radial component of the center-of-mass velocity and we have used the fact that center of mass only experiences a radial acceleration due to the third body. Thus we get \[\delta \psi_\mathrm{G} \approx -2\frac{Gm_3}{r^2c} \frac{V_r}{c} \frac{t^2}{2} (2\pi f) \: .\] The ratio \(\mathcal{R}_{\rm DG}\) between gravitational redshift and Doppler induced phase shifts : \[\begin{align} \lvert\mathcal{R}_{\rm DG}\rvert \approx 2\frac{V_{\rm r}}{c P}, \end{align}\] where we neglected a factor of order unity. We see that the Doppler shift due to gravitational redshift is formally suppressed compared to the Doppler shift due to line-of-sight motion by a factor \(V_{\rm r}/c\). The magnitude of the radial velocity of binary on an eccentric orbit around the third body scales approximately as: \[\begin{align} V_{\rm r} \sim c \left(\frac{r_{\rm S}}{r_0} \right)^{1/2}e_{3}, \end{align}\] where \(e_3\) is the orbital eccentricity and \(r_{\rm S}\) is the Schwarzschild radius of the third body. For binaries around massive black holes, such radial velocities could therefore easily reach order \(0.01c\). Therefore, while we expect Doppler dephasing to dominate, dephasing due to radial motion in the third body potential may dominate for small values of the projection factor.

Before moving on to eccentric waveforms let us note that analytic formulas for the magnitude of Doppler-dephasing in eccentric binaries have already been worked out in some recent studies [49]–[51] and further elaborated on by Ref. [36].

5 Dephasing in eccentric waveforms↩︎

On a conceptual level, the presence of orbital eccentricity adds several new aspects to the analyisis of dephasing perfomred in previous sections. Firstly, eccentricity influences the functional dependence of the GW phase, and therefore of any dephasing, on the GW frequency. Additionally, other ways in which perturbations can modify the orbital evolution emerge, related to the extra degrees of freedom the binary orbit has. We list here a few considerations:

• An external angular momentum flux can cause a secular perturbation of the binary eccentricity in addition to affecting the semi-major axis and the orbital frequency itself. This can indirectly induce a modification to the chirp of the signal.

• Any precession of the argument of pericentre induced by a perturbing force directly influences the GW phase and modulates the amplitudes of the GW polarisations. Additionally, changes in the binary orientation due to precession will also in turn change the coupling between the binary and the perturbing force. In fact, 1PN forces will already necessarily make the orbit precess, which precession can then interact with perturbing forces in non-trivial ways. Furthermore, the misalignment of spin angular momentum of the binary components with the orbital angular momentum of the binary leads to the precession of the orbital plane, further modulating the GW signal and making off-plane forces interact with the orbit.

• The presence of external perturbative potientials breaks the Keplerian definitions of both semi-major axis and eccentricity in non trivial ways, similarly to the discussion in section 3.3.

Here we do not aim to cover all of these options in full. Instead, we focus on how a generic dephasing prescription can be translated between time and frequency domains, analogously to the case for cricular orbits in section 3.1. Additionally, we analyse the behaviour of dephasing prescriptions that arise from energy and angular momentum fluxes, i.e. a generalised version of section 3.2.

5.1 Basics of eccentric waveforms↩︎

We start our discussion with a review of the basics of eccentric waveforms, based on the seminal work of Yunes et al. [116]. While at a given time instance circular binaries only emit quadrupolar GWs at a singular frequency, i.e. twice the orbital frequency, eccentric binaries have a wide spectrum corresponding to the various epicyclic frequencies of the orbit, i.e. multiples of the orbital frequency. This makes the GW strain at a given time a sum of of different harmonics: \[h(t) = \sum_{\ell=1}^{\infty} \mathcal{A}_\ell(t) \exp\left[ -i\ell \Phi(t) \right] \: , \label{eq:Ecc95waveform}\tag{19}\] where \(\Phi\) is the mean anomaly, while \(\ell\) corresponds to the epicyclic expansion. Fig 4 illustrates this expansion for a single period of an orbit with \(e=0.5\). Note that \(\mathcal{A}_\ell\) now is slowly varying complex quantity, the phase of which depends on the projection of the elliptic orbit on the observer’s frame and thus can vary due to the precession of the orbit. Such an expansion for non-precessing binaries was carried out for low eccentricities by [116]. The mean anomaly in Eq. 19 will evolve as \[\Phi(t) = \int^t 2\pi F(t') \, \mathrm{d}t' \: ,\] where \(F\equiv \Omega/2\pi\) is the orbital frequency. Then in a similar manner as we did in the circular case we can use the stationary phase approximation to translate Eq. 19 into a Fourier domain waveform. This can be done as long as the two conditions for the approximation are valid, i.e. \(\mathrm{d}(\ln \mathcal{A}_\ell)/\mathrm{d}t \ll \ell\mathrm{d}\Phi/\mathrm{d}t\), while the frequency is also evolving slowly, i.e. \(\ell\mathrm{d}^2\Phi/\mathrm{d}t^2 \ll \ell^2(\mathrm{d}\Phi/\mathrm{d}t)^2\), where we replaced \(\phi\) by \(\phi_\ell = \ell \Phi\). Note that in general, the phase of the Fourier integrand can encounter catastrophes, i.e locations where multiple stationary points coalesce and the second time derivative of the phase of the integrand vanishes together with the first derivative [117], [118]. Such catastrophes render the stationary phase approximation ill-suited for the evaluation of the Fourier integral, and modifications to this approximation scheme must be employed [119]–[121]. In practice, as long as most of the power of the eccentric waveform comes from harmonics below \(\ell_\mathrm{max}\), where catastrophes are not present, the stationary phase approximation remains valid. Then we can write:

\[\tilde{h}(f) \approx \sum_{\ell=1}^{\ell_\mathrm{max}} \tilde{\mathcal{A}}_\ell\left( f \right) \exp\left[ i\psi_\ell(f) \right] = \sum_{\ell=1}^{\ell_\mathrm{max}} \tilde{h}_\ell(f) \: ,\] where \(f=\ell F\) is still the GW frequency. The phases \(\psi_\ell(f)\) are just multiples of the Fourier phase: \[\psi_\ell (f) = \ell \Psi\left( F=\frac{f}{\ell} \right) \: ,\] where \[\Psi(F) = 2\pi F t(F) - \Phi(t_F) \: .\] The formulas for \(t(F)\) and \(\Phi(t_F)\) are then also analogous to the circular case: \[\begin{align} t(F) &= \int^F \dot{F}^{-1} \, \mathrm{d}F' \: , \tag{20} \\ \Phi(t_F) &= 2\pi \int^F F' \dot{F}^{-1} \, \mathrm{d}F' \tag{21} \: . \end{align}\] Hence, there is a one-to-one mapping between \(\Psi\) and the orbital frequency \(F\), independent of \(\ell\), which is not the case for \(\Psi\) and \(f\). The Fourier phase \(\Psi\) at a specific \(F\) will give contribution to phases at multiple different GW frequencies \(f=\ell F\), or equivalently, the Fourier domain waveform and thus its amplitude and phase at \(f\) will be accumulated from different stages of the binary evolution for different \(\ell\) that correspond to an orbital frequency \(F/\ell\). The Fourier domain amplitude of the different harmonics and the total waveform is illustrated in Fig. 5. While the amplitude of each of the harmonics is smooth, note how when added up they create a noisy waveform due to their incoherent phases.

Also in an analogous way to the circular case we have the following relation between Fourier phase and frequency evolution: \[\begin{align} \frac{\mathrm{d}\Psi}{\mathrm{d}F} &= 2\pi t(F)\: , \tag{22}\\ \frac{\mathrm{d}^2\Psi}{\mathrm{d}F^2} &= \frac{2\pi}{\dot{F}}\: . \tag{23} \end{align}\]

While so far the formulas have been analogous to the circular case, for eccentric binaries the frequency evolution is modified in a non-trivial way: \[\begin{align} \dot{F} = \dot{F}^{e=0}\times \mathcal{F}(e) \end{align}\] where \[\dot{F}^{e=0}(F) = \frac{96}{5}(2\pi)^{8/3} \mathcal{M}^{5/3}F^{11/3} \: , \label{eq:dFvace0}\tag{24}\] and the enhancement function \(\mathcal{F}(e)\) reads [122]: \[\begin{align} \mathcal{F}(e) = \left(1+\frac{73}{24}e^2 + \frac{37}{96}e^4 \right)(1-e^2)^{-7/2} \: . \label{eq:Fe} \end{align}\tag{25}\] In order to evaluate the integrals in Eqs. 20 and 21 we need to have the relation \(e(F)\). While \(F(e)\) can be given as an analytic function, the inverse formula generally cannot be given in analytic form. However, for small eccentricities, a series expansion can be carried out. Since the aim of this work is only to illustrate different types of dephasings, we only carry out the expansion to leading order in eccentricity, however, expansions of up to \(\mathcal{O}(e^9)\) can be found for example in Ref. [116]. The relation between the semi-major axis and the eccentricity with the 2.5PN GW radiation term reads: \[a(e) = \left( \frac{M}{4\pi^2F^2} \right)^{1/3} = \frac{a_\mathrm{in}}{g(e_\mathrm{in})}g(e) \: , \label{eq:ae}\tag{26}\] where \(e_\mathrm{in}\) is the initial eccentricity at the initial semi-major axis \(a_\mathrm{in}\), and \[\begin{align} g(e) = \frac{e^{12/19}}{1 - e^2}\bigg(1 + \frac{121}{304}e^2\bigg)^{870/2299} \: . \label{eq:ge} \end{align}\tag{27}\] Then to leading order \(e(F)\) reads: \[e(F) \approx e_\mathrm{in} \left( \frac{F}{F_\mathrm{in}} \right)^{-19/18} \: ,\] where \(F_\mathrm{in}\) is the initial orbital frequency.

5.2 Time-domain and Fourier-domain dephasing↩︎

After finding an appropriate representation of the GW signal in the Fourier domain our task then reduces to finding \(\delta\Psi(F)\), from which the dephasing for all the different harmonics can be found by \(\delta\psi_\ell(f/\ell) = \ell\delta\Psi(F)\) (see illustration in Fig. 6). Most of the statements that we found for the GW dephasing in the circular case remains valid for \(\delta\Psi(F)\) as well. This means that assuming a perturbation in the chirp of the orbital frequency: \[\dot{F}_\mathrm{tot}(F) = \dot{F}_\mathrm{vac}(F) + \delta \dot{F}(F) \: ,\] we get that \[\begin{align} \delta t(F) &\approx- \int^F\frac{\delta\dot{F}}{{\dot{F}_\mathrm{vac}}^2} \, \mathrm{d}F' \: , \\ \delta\varphi(F) &= \Phi(t_F) - \Phi_\mathrm{vac}(t_{F,\mathrm{vac}}) \nonumber \\ &\approx - 2\pi \int^F\frac{F'\delta\dot{F}}{{\dot{F}_\mathrm{vac}}^2} \, \mathrm{d}F' \: . \end{align}\] And thus: \[\delta\Psi(F) = 2\pi F\delta t(F) - \delta\varphi(F) \: .\] For the second derivative with respect to frequency we get: \[\frac{\mathrm{d}^2(\delta\Psi)}{\mathrm{d}F^2} = - 2\pi \frac{\delta \dot{F}}{\dot{F}_\mathrm{vac}^2} \: .\] Note that in this case as well \(\delta\Psi(F)\) itself is free up to two integration constants: \[\delta\Psi \rightarrow \delta\Psi + 2\pi F t_c - \Phi_c \: .\] The relation between the time- and Fourier domain dephasings also holds: \[\delta \Phi(t_F) \approx -\delta\Psi(F) \: ,\] meaning that in order to get the Fourier dephasing at an orbital frequency \(F\), to first order it is sufficient to evaluate the difference in the mean anomalies of the perturbed and unperturbed orbits at \(t_F\) when the orbital frequency of the unperturbed orbit reaches \(F\). This relation becomes even more useful in the eccentric case, as the analytic evaluation of the exact dephasing can become complicated, and so direct numerical simulations can be used instead to obtain \(\delta\Phi(t_F)\). Then \(\delta\Psi(F)\) can be used to calculate \(\delta\)SNR with an analytic waveform, or simply to validate analytic dephasing prescriptions.

All this means that the formulas in columns \(3-5\) of Table 1 for intrinsic and Doppler dephasings remain valid for eccentric orbits as well by replacing \(f\), \(\phi\) and \(\psi\) with \(F\), \(\Phi\) and \(\Psi\), respectively. The exact form of \(\delta\dot{F}(F)\) can however be complicated, as we will see later, due to the mutual dependence of the parameters on each other.

5.3 Dephasing due to additional energy and angular momentum fluxes↩︎

An additional energy flux generated by an external perturbation has a similar direct effect as in the circular case: \[\delta \dot{F} \approx \dot{F}_\mathrm{vac} \frac{\delta P}{\dot{E}_\mathrm{vac}} \: .\] However, in the case of eccentric orbits we have an additional, indirect effect as well. This arises from the fact that the chirp is not only a function of the semi-major axis (or equivalently, the orbital frequency) but the eccentricity as well. Thus if due to an external perturbation the eccentricity at a given orbital frequency \(F\) is modified with respect to what the leading order evolution from the 2.5PN radiation would dictate, the chirp will also be modified as: \[\delta \dot{F} \approx \frac{\partial \dot{F}_{\rm{vac}} }{\partial e}\delta e \: ,\] where \[\begin{align} \frac{\partial \dot{F} }{\partial e} &= \dot{F}^{e=0}\left(\frac{307}{23}e+\frac{9611}{552}e^3+\frac{37}{32}e^5 \right)\left(1-e^2 \right)^{-9/2} \end{align}\] Figure 7 illustrates the cause of this extra dephasing. The solid lines show different evolutionary paths the binary can take depending on the initial eccentricity \(e_\mathrm{in}\) at some initial semi-major axis \(a_\mathrm{in}\). A perturbation in eccentricity or in semi-major axis can divert the binary from its original trajectory, and put it on a trajectory corresponding to an initial eccentricity \(e_\mathrm{in}+\delta e_\mathrm{in}\). Thus the eccentricity at a given orbital frequency will also be modified by \(\delta e(F)\). A surprising result is that even though a change in the semi-major axis will not explicitly modify the eccentricity, by changing the orbital frequency while keeping the eccentricity constant it will put the binary on a different trajectory and thus the eccentricity after the frequency change will differ from the expected eccentricity. A perturbation will only keep the binary on the same trajectory if \(\delta\dot{e}/\delta\dot{a}\propto \dot{e}_\mathrm{vac}/\dot{a}_\mathrm{vac}\).

Parameterizing the trajectories with \(e_\mathrm{in}\) and differentiating Eq. 26 with respect to time we have that for general \(\dot{a}\) and \(\dot{e}\): \[\dot{a} = -\frac{a_\mathrm{in}}{g(e_\mathrm{in})^2}g'(e_\mathrm{in}) g(e) \dot{e}_\mathrm{in} + \frac{a_\mathrm{in}}{g(e_\mathrm{in})} g'(e) \dot{e} \: ,\] where \(g'(e)\equiv \mathrm{d}g(e)/\mathrm{d}e\). Note that the 2.5PN radiation terms give \(\dot{e}_\mathrm{in}=0\) and thus from the perturbations we get: \[\dot{e}_\mathrm{in} = \frac{g(e_\mathrm{in})}{g(e)} \frac{g'(e)}{g'(e_\mathrm{in})} \delta \dot{e} - \frac{g(e_\mathrm{in})^2}{g(e)g'(e_\mathrm{in})} \frac{\delta\dot{a}}{a_\mathrm{in}} \: .\] Then the total change in \(e_\mathrm{in}\) parameterizing the trajectory will be: \[\delta e_\mathrm{in}(F) = \int^{t_F} \dot{e}_\mathrm{in} \, \mathrm{d}t' \: .\] We then need to determine what this change in the initial eccentricity parameter will mean for \(\delta e(F)\). Differentiating the relation \(a(e,e_\mathrm{in})\) for fixed semi-major axis to linear order we get that a change \(\delta e_\mathrm{in}\) results in \[\delta e \approx \frac{g(e)}{g(e_\mathrm{in})} \frac{g'(e_\mathrm{in})}{g'(e)} \delta e_\mathrm{in} \: .\] Combining these relations we can now determine how perturbations in eccentricity and semi-major axis modify the eccentricity at a given orbital frequency. Replacing the time integral with a frequency integral and keeping everything to linear order we have: \[\begin{align} \delta e(F) = \frac{g(e_0)}{g'(e_0)}\int^F \frac{{\rm d}F'}{\dot{F}_{\rm vac}} \left[ \frac{g'(e_0)}{g(e_0)} \delta\dot{e} - \frac{\delta\dot{a}}{a_0}\right] \: , \end{align}\] where \(e_0 \equiv e_\mathrm{vac}(F)\) and \(a_0 \equiv a_\mathrm{vac}(F)\) and \(F'\) is a dummy integration variable. At lowest order in eccentricity this reduces to: \[\begin{align} \delta e(F) = e_0(F)\int^F \frac{\delta\dot{e}(F')}{\dot{F}_{\rm vac}^{{e=0}}(F') e_0(F')} \, {\rm d}F' \: . \end{align}\] Thus, perturbations in the semi-major orbit do not modify the trajectory up to linear order in eccentricity. This can be understood by the fact that the trajectories in Fig. 7 are almost horizontal for small eccentricities, thus changing the the semi-major axis would keep you on the same trajectory to leading order. Along the same argument, a small change in the eccentricity changes the initial eccentricity by an increasingly large amount as the eccentricity decreases.

These eccentricity and semi-major axis perturbations can be induced by external energy and angular momentum fluxes (torques, \(\delta\tau\)). These quantities are related by the following relations: \[\begin{align} \delta \dot{e} &= \frac{e^2-1}{e} \left( \frac{\delta\tau}{L_\mathrm{orb}} + \frac{\delta P}{2 E_\mathrm{orb}} \right) \: , \\ \delta \dot{a} &= -a \frac{\delta P}{E_\mathrm{orb}} \: , \end{align}\] where \(\delta P\) and \(\delta \tau\) are negative when the binary is losing energy or angular momentum.

Thus the induced eccentricity finally becomes \[\begin{align} \delta e(F) = \frac{g(e_0)}{g'(e_0)}&\int^F \frac{{\rm d}F'}{\dot{F}_{\rm vac}} \left[ \frac{g'(e_0)}{g(e_0)} \frac{e_0^2-1}{e_0} \frac{\delta \tau}{L_\mathrm{orb}} \right. \nonumber \\ &\left. + \left(\frac{g'(e_0)}{g(e_0)} \frac{e_0^2-1}{2e_0}+1\right) \frac{\delta P}{E_\mathrm{orb}} \right] \: . \end{align}\] To lowest order in eccentricity this means \[\begin{align} &\delta e(F) =\nonumber\\& -e_0(F)\int^F \frac{1}{\dot{F}_{\rm vac}^{{e=0}} e_0^2} \left( \frac{\delta\tau}{L_\mathrm{orb}} + \frac{\delta P}{2 E_\mathrm{orb}} \right) \, {\rm d}F' \: , \end{align}\] and once again the total induced perturbation to the frequency chirp is²: \[\begin{align} \delta \dot{F}_{\rm tot} = \dot{F}_{\rm vac}\frac{\delta P}{\dot{E}_{\rm vac}} + \frac{\partial \dot{F}_{\rm vac}}{\partial e}\delta e. \label{eq:Fd95ecc95twoterms} \end{align}\tag{28}\] These results of our calculations are summarized in Table 2.

It is interesting to ask which of the two components of the dephasing dominates. While of course the answer depends on the specifics of the perturbing force, we work through an illustrative example where: \[\begin{align} \frac{\delta\dot{a}}{a} &= {\rm const.} \equiv K_{\rm a}\\ \delta\dot{e} &=0. \end{align}\] where one would naively expect no dephasing contribution to dephasing from an eccentricity perturbation. We can calculate the resulting dephasing via the usual prescriptions reported in Table 2. here we report the time-domain binary dephasing at a given frequency, which is given by the integral: \[\begin{align} \delta \varphi(F) = 2 \pi\int {\rm d}F' \left[F'\frac{\delta \dot{F}_{\rm tot}}{\dot{F}_{\rm vac}^2} \right]\equiv \delta \varphi_a + \delta \varphi_e \end{align}\] where \(\delta \varphi_a\) comes from the first term in Eq. 28 and \(\delta \varphi_e\) from the second. We can solve the integral analytically in the limit of small eccentricity, finding that the ration of the two induced phase-shifts reads: \[\begin{align} \frac{\delta \varphi_e}{\delta \varphi_a} \approx - \frac{116337}{14848}e_0(F)^2 \approx- 7.8\times e_0(F)^2. \end{align}\] Therefore, even without any specific eccentricity perturbation \(\delta \dot{e}\), the dephasing due to the induced change in eccentricity \(\delta e\) can still be dominant for \(e_0 \sim 0.4\). Compare now a similar calculation in which \(\delta\dot{e} = K_{e}\) is constant. Repeating the integral yields: \[\begin{align} \frac{\delta \varphi_e}{\delta \varphi_a} \approx - \frac{55107}{5626} \frac{K_{\rm e}}{K_{\rm a}}e_0(F). \end{align}\] That is, the dephasing due to the induced change in eccentricity can already dominate for \(e_0\sim 0.1\), provided that \(\delta \dot{a}/a\) and \(\delta \dot{e}\) are of similar order.

Another interesting and more physically relevant example would be that of a von Zeipel–Kozai–Lidov (ZKL) eccentricity oscillation induced by a tertiary object orbiting around the binary [123]–[125]. In such a case in the secular limit the inner binary exchanges angular momentum without exchanging energy, thus an eccentricity perturbation is induced without perturbing the inner binary’s semi-major axis. For a circular outer binary and optimal relative orientation of the inner and outer binaries we have that (see e.g. [126]): \[\frac{\delta \dot{e}}{e} \approx \frac{15}{8} \sqrt{1-e^2} \frac{T_\mathrm{in}}{T_\mathrm{out}^2} \: ,\] where \(T_\mathrm{in}\) is the orbital period of the inner binary and \(T_\mathrm{out}\) is that of the outer binary formed by the tertiary object and the inner binary. In the regime where ZKL oscillations dominate over GW radiation, it has been shown that eccentricity oscillations can cause a gradual phase shift over many oscillation periods, which might be detectable by LISA. In our case, where we are considering perturbations that are small compared to the 2.5PN radiation reaction terms, 1PN precession of the binary will also inevitably become important, which will suppress ZKL oscillations. However, for pedagogical reasons let us now imagine an environmental effect that causes an eccentricity perturbation with the same scaling with eccentricity as the ZKL effect. Then by integrating Eq. 28 once again we can find that the scaling of the resulting induced dephasing would be: \[\delta \varphi_{\rm e} \propto e_{0}(F)^2 F^{-16/3}.\]

6 SNR and \(\delta\)SNR calculations↩︎

We conclude our overview of dephasing with a brief discussion of signal-to-noise (SNR) ratios, which are related to the observability of a GW signal and any potential perturbation. The SNR is calculated from the Fourier domain strain according to: \[\mathrm{SNR}^2 = 4\int \frac{|\tilde{h}(f)|^2}{S_n(f)} \mathrm{d}f \: ,\] where \(S_n(f)\) is the spectral noise density curve of the detector. The interpretation of this integral is that of summing the excess power of the GW signal with respect to the detector noise, in every frequency bin. Similarly, the observability of a perturbation is usually estimated by the SNR of the difference signal [6], i.e. \[\delta \mathrm{SNR}^2 = 4 \int \frac{|\delta\tilde{h}(f)|^2}{S_n(f)} \mathrm{d}f \: ,\] where \[\begin{align} \delta \tilde{h}(f) &= \tilde{h}_\mathrm{tot}(f) - \tilde{h}_\mathrm{vac}(f) \approx \tilde{h}_\mathrm{vac}(f)\left( e^{i\delta\psi(f)}-1 \right) \nonumber\\ &\approx \tilde{h}_\mathrm{vac}(f)\cdot i\delta\psi(f) \: , \end{align}\] where we assumed that the main contribution to the difference signal comes from the dephasing rather than the amplitude modulation. Thus \[\delta \mathrm{SNR}^2 \approx 4 \int |\delta\psi(f)|^2 \frac{|\tilde{h}(f)|^2}{S_n(f)} \mathrm{d}f \: .\]

The SNR calculations for eccentric sources is also analogous to the circular case, only in this case the total SNR comes from a contribution from all the different harmonics: \[\mathrm{SNR}^2 = 4\int \frac{1}{S_n(f)} \Big|\sum_{\ell}\tilde{h}_\ell(f) \Big|^2 \mathrm{d}f \: . \label{eq:SNR95ecc}\tag{29}\] Then the difference signal is: \[\begin{align} \delta \tilde{h}(f) &\approx \sum_\ell \tilde{h}_\mathrm{\ell,vac}(f)\cdot i\delta\psi_\ell(f) \: , \label{eq:deltah95ecc} \end{align}\tag{30}\] and thus the \(\delta\)SNR: \[\delta \mathrm{SNR}^2 \approx 4\int \frac{1}{S_n(f)} \Big|\sum_{\ell} \delta \psi_\ell(f)\tilde{h}_\ell(f) \Big|^2 \mathrm{d}f \: ,\] where \[\delta \psi_\ell(f) = \ell \delta\Psi(f/\ell) \: .\] Note that the dephasing for harmonics with \(\ell\gg 1\) can easily become large and even exceed \(2\pi\), in which case the approximation in Eq. 30 clearly does not hold. Therefore the use of this formula requires some consideration and \(\exp[i\delta\psi_\ell(f)]-1\) should be used instead of \(i\delta\psi_\ell(f)\) in the general case.

Due to the incoherence of the phases of different harmonics at a given GW frequency \(f\) (see Fig 4), we can also assume that different harmonics give orthogonal contributions to the total SNR [127], i.e. the cross terms in the sum in Eq. 29 can be neglected and thus the sum can be moved outside of the absolute value operator: \[\mathrm{SNR}^2 \approx 4\int \sum_{\ell}\frac{|\tilde{h}_\ell(f) |^2}{S_n(f)} \mathrm{d}f \: .\] We can also make a similar assumption for the orthogonality of different harmonics that contribute to \(\delta\)SNR: \[\delta\mathrm{SNR}^2 \approx 4\int \sum_{\ell}\frac{ |\delta\tilde{h}_\ell(f) |^2}{S_n(f)} \mathrm{d}f \: .\]

We note that defining the detectability in terms of the \(\delta\)SNR can be useful to explore large parameter spaces efficiently, as the generation of accurate waveforms can be expensive and \(\delta\)SNR values instead provide fast detectability criteria. However, this way of defining detectability does not take into account possible degeneracies between parameters and can be too generous in estimating the detectability of an effect. Although potentially time-consuming, other, more accurate methods exist for estimating detectability or waveform accuracy, such as the recent update on the mismatch criterion in Ref. [128].

6.1 \(\delta\)SNR enhancement in eccentric binaries↩︎

Higher harmonics of the GW signal can contribute significantly to the total SNR of eccentric binaries, which can make the SNR of eccentric binaries higher in some mass and initial eccentricity range [129]. This enhancement with eccentricity can be even more pronounced for the \(\delta\)SNR of dephased eccentric binaries. To understand this, let us consider the contribution of the \(\ell=2\) mode, which is dominant for circular binaries, and that of the \(\ell=4\) mode, taken at a given GW frequency \(f\). Let us suppose that for \(f_{\ell=2} = 2F = f\), the phase shift is \(\delta\psi_{\ell=2}(f)=2\delta\Psi(f/2)\equiv\delta\psi\). The contribution for the same GW frequency \(f\) from the \(\ell=4\) comes from an earlier stage of the binary evolution. As an example, let us now consider a simple power-law scaling of \(\Psi(F)\propto F^{-n}\), where \(n=13/3\) corresponds to phase shifts due to Roemer-delay in the circular case. With this scaling, the phase shift of the \(\ell=4\) at \(f=4F\) is \(\delta \psi_{\ell=4}(f) = 4\delta\Psi(f/4) = 2\cdot 2^{n} \delta\psi\). Thus we get a factor of \(\sim40\) enhancement for the contribution to the \(\delta\)SNR at GW frequency \(f\) from the \(\ell=4\) mode. This huge enhancement factor can make higher harmonics dominate the \(\delta\)SNR even for low eccentricities, where the SNR contribution of these harmonics is normally suppressed (see Fig. 5.

Let us now be more exact and calculate the \(\delta\)SNR from Roemer-delay for eccentric binaries. The phase shift due to line-of-sight acceleration of the center of mass of an inspiralling binary results in a phase shift [50]: \[\begin{align} \left|\delta\Psi(F)\right| &\approx \left|\delta\Phi(t_F)\right| \approx \frac{2\pi\Delta t}{T} \nonumber \\ &= 2\pi \frac{a_\mathrm{CM}}{2c} t_F^2 F \: , \end{align}\] where \(a_\mathrm{CM}\) is the projection of the center-of-mass acceleration in the line-of-sight direction. The inspiral time for circular binaries, as it was shown in Sec. 2: \[t_F^\mathrm{circ} = 5(16\pi F)^{-8/3} \mathcal{M}^{-5/3} \: .\] For non-zero eccentricities the inspiral time from a given orbital frequency is reduced. According to Ref. [130] it can be estimated by the following formula: \[t_F^e \approx t_F^\mathrm{circ} \frac{R(e)}{\mathcal{F}(e)} \: ,\] where \(\mathcal{F}(e)\) is defined by Eq. 25 and \(R(e) = 8^{1-\sqrt{1-e}}\). Even though the extra factor reduces the inspiral time significantly, it is still \(\sim0.8\) for \(e=0.3\). Therefore, for eccentricities \(e\lesssim0.3\) the eccentric inspiral time is approximately the circular inspiral time. Thus, for simplicity we assume a dephasing prescription of \(\delta\Psi(F) = C \cdot F^{n}\) with \(n=-13/3\), and investigate the change in the \(\delta\)SNR with increasing orbital eccentricity at some reference frequency. We leave the detailed investigation of \(\delta\)SNR for higher eccentricities to follow-up studies. We choose the reference frequency as the frequency of the \(\ell=2\) mode at the lower cut-off frequency of a given GW detector. For the frequency-domain waveform, we use the model of Ref. [131] without 1PN precession.

The results are shown in Fig. 8. We can see that even though the total SNR of binaries has a weak dependence on the eccentricity for moderate eccentricities, the \(\delta\)SNR from the phase shift increases rapidly with eccentricity and obtains a factor of \(\sim5\) enhancement for \(e_\mathrm{10Hz}\approx0.4\) for LIGO and a factor of \(\sim10^3\) enhancement for ET at \(e\approx0.3\). Interestingly but not surprisingly, apart from close to zero eccentricities, the contributions to the \(\delta\)SNR is dominated by harmonics \(\ell>2\). This is even more pronounced for the results for ET, where already at \(e_\mathrm{2Hz}\gtrsim0.2\) the contributions from \(\ell>10\) harmonics dominate. Note that this is not the case for the total SNR, where, except for extreme eccentricities, the \(\ell=2\) mode always dominates. Note that for \(\ell\gg1\), large contributions to \(\delta\)SNR can come from early stages of the binary evolution, where \(e\gtrsim0.9\), and thus our power-law phase shift prescription is no longer valid. Hence, our calculations here represent an idealistic scenario. For large eccentricities, the \(\delta\)SNR is expected to be reduced, however, the significant amplification at \(e\sim0.1\) will remain valid. Note also that calculating the \(\delta\)SNR this way does not take into account possible degeneracies with other parameters and also assumes that these orbital parameters are known exactly. For real observations this will obviously not be the case. Nevertheless, these results underline the significance of eccentric binaries in the search for environmental effects in GW signals.

7 Conclusion↩︎

Here we provide a short summary of the main conclusions of this paper and emphasize the main take-aways. In Section 2, we went through the derivation of the Fourier domain waveform using the stationary phase approximation for circular binary inspirals, based on the seminal work of Cutler & Flanagan. In particular we highlighted how the phases in the time and Fourier domains are related. Section 3 presented the relevant expressions and highlighted the key differences between various types of dephasing across the two domains. We also outlined the primary environmental effects that cause phase shifts in circular binaries, including those arising from extra energy losses, radial forces, and Doppler motion. The most important formulas are reported in Table 1. In Section 4, we explored typical astrophysical scenarios where multiple sources of dephasing occur simultaneously. We assessed their relative impacts with simple but informative toy models. In Section 5 we extended this analysis to eccentric binaries. The presence of eccentricity introduces a variety of ways in which environmental perturbations can induce dephasing. We focused here on the ways in which additional energy and angular momentum fluxes can influence the binary inspiral, deriving a novel prescription for the gravitational wave phase. Finally, in Section 6, we briefly addressed how signal-to-noise ratios affect the ability to detect these environmental dephasing effects, and showed in an idealized example how the SNR from the difference signal due to environmental phase shifts can be significantly amplified for non-zero eccentricities.

To conclude, we delineate again the aim of this paper is twofold. Firstly, it serves as a cohesive overview of the topic of GW dephasing, clarifying the underlying concepts and organizing essential formulas, examples, and methodological tools in a practical format. It aims to support ongoing research while also guiding newcomers through a modern and accessible entry into the field. Secondly, it highlights the interesting possibilities of analyzing dephasing prescriptions in eccentric GW signals, as the foundations worked out here already show how interesting features arise from the interplay between semi-major axis and eccentricity evolution.

References↩︎