<figure><img src="image_1.png" /></figure>

# Gravitational-wave Detection and Parameter Estimation for Accreting Black-hole Binaries and Their Electromagnetic Counterpart

Andrea Caputo¹, Laura Sberna², Alexandre Toubiana³,⁴, Stanislav Babak³,⁵, Enrico Barausse⁴,⁶,⁷,⁸, Sylvain Marsat³, and Paolo Pani⁹

¹Institut de Física Corpuscular—CSIC/Universitat de València, Parc Científic de Paterna. C/Catedrático José Beltrán, 2 E-46980 Paterna (Valencia), Spain

²Perimeter Institute, 31 Caroline St N, Ontario, Canada

³APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité, 10, rue Alice Domon et Léonie Duquet F-75205 Paris Cedex 13, France

⁴Institut d'Astrophysique de Paris, CNRS & Sorbonne Universités, UMR 7095, 98 bis bd Arago, F-75014 Paris, France

⁵Moscow Institute of Physics and Technology, Dolgoprudny, Moscow region, Russia

⁶SISSA, Via Bonomea 265, I-34136 Trieste, Italy

⁷INFN, Sezione di Trieste, Italy

⁸IFPU—Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

⁹Dipartimento di Fisica, "Sapienza" Università di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 5, I-00185, Roma, Italy

Received 2020 January 18; revised 2020 February 27; accepted 2020 February 27; published 2020 April 1

## Abstract

We study the impact of gas accretion on the orbital evolution of black-hole binaries initially at large separation in the band of the planned Laser Interferometer Space Antenna (LISA). We focus on two sources: (i) stellar-origin black-hole binaries (SOBHBs) that can migrate from the LISA band to the band of ground-based gravitational-wave (GW) observatories within weeks/months; and (ii) intermediate-mass black-hole binaries (IMBHBs) in the LISA band only. Because of the large number of observable GW cycles, the phase evolution of these systems needs to be modeled to great accuracy to avoid biasing the estimation of the source parameters. Accretion affects the GW phase at negative (−4) post-Newtonian order, being thus dominant for binaries at large separations. Accretion at the Eddington or at super-Eddington rate will leave a detectable imprint on the dynamics of SOBHBs. For super-Eddington rates and a 10 yr mission, a multiwavelength strategy with LISA and a ground-based interferometer can detect about 10 (a few) SOBHB events for which the accretion rate can be measured at 50% (10%) level. In all cases, the sky position can be identified within much less than 0.4 deg² uncertainty. Likewise, accretion at $\gtrsim$100% of the Eddington rate can be measured in IMBHBs up to redshift $z \approx 0.1$, and the position of these sources can be identified within less than 0.01 deg² uncertainty. Altogether, a detection of SOBHBs or IMBHBs would allow for targeted searches of electromagnetic counterparts to black-hole mergers in gas-rich environments with future X-ray detectors (such as Athena) and/or radio observatories (such as SKA).

*Unified Astronomy Thesaurus concepts: Astrophysical black holes (98); Gravitational waves (678); Accretion (14)*

## 1. Introduction

Among the main gravitational-wave (GW) sources detectable by the future Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) are binary black holes with relatively small masses, down to a few tens of solar masses (Sesana 2016). LISA can detect these systems when they are still at large separations and thus probe their low-frequency dynamics. In more detail, these systems include: (i) stellar-origin black-hole binaries (SOBHBs) of a few tens up to $\sim 100 M_{\odot}$, whose coalescences are also observed by terrestrial GW detectors (Abbott et al. 2019); and, if they exist, (ii) intermediate-mass black-hole binaries (IMBHBs) with component masses in the range ($10^2, 10^5) M_{\odot}$ (Miller & Colbert 2004).

SOBHBs will be first observed in the LISA $\sim$ mHz band, and will then disappear for weeks/months before entering the $\gtrsim$1 Hz band of ground detectors, where they merge (Sesana 2016). Despite this frequency gap, piecing together the LISA low-frequency regime and the terrestrial high-frequency merger will allow for effectively observing these systems for $10^5-10^6$ GW cycles. Therefore, even small inaccuracies in modeling the GW phase evolution will bias the estimation of the parameters (and particularly the merger time) or even prevent detection by LISA.

IMBHBs might be detected by LISA for the first time for a whole range of total masses and mass ratios, with the lighter binaries spending more time in band. While the existence of

intermediate-mass black holes has not been confirmed yet, several candidates exist (see, e.g., Mezcua 2017, for a review), and they might also provide seeds for the growth of the supermassive black holes that are ubiquitously observed in the local universe (see, e.g., Latif & Ferrara 2016; Mezcua 2017). While their formation mechanism is unknown, proposed scenarios include direct collapse of massive first-generation, low-metallicity Population III stars (Madau & Rees 2001; Schneider et al. 2002; Kinugawa et al. 2014; Ryu et al. 2016), runaway mergers of massive main-sequence stars in dense stellar clusters (Miller & Hamilton 2002; Portegies et al. 2002, 2004; Atakan Gurkan et al. 2004; Mapelli 2016), accretion of residual gas onto stellar-origin black holes (Leigh et al. 2013), and chemically homogeneous evolution (Marchant et al. 2016).

Both SOBHBs and IMBHBs offer the potential to constrain low-frequency modifications of the phase evolution, if the latter are included in the GW templates used for the analysis in the LISA band. Such low-frequency phase modifications may appear, e.g., if the dynamics of these systems is governed by a theory extending/modifying general relativity (Barausse et al. 2016; Carson & Yagi 2020a; Gnocchi et al. 2019), or as a result of interactions (already within general relativity) of the binary with the surrounding gas, if the latter is present (Barausse et al. 2014, 2015; Cardoso & Maselli 2019; Tamanini et al. 2020).
There is currently no evidence that the SOBHBs observed by GW detectors live in gas-rich environments—and no electromagnetic (EM) counterpart to these sources has been detected so far (Abbott et al. 2016b). Binaries involving accreting stellar-origin black holes are observed in X-rays (Charles & Coe 2003), but the accreting gas is provided by a stellar companion. However, gas may be present earlier in the evolution of SOBHBs, and some of it may survive in the binary’s surroundings. For instance, in the field-binary formation scenario (Abbott et al. 2016a) for SOBHBs, gas plays a key role in the common envelope phase, although the latter typically precedes the merger by several Myr. Also note that SOBHBs may form preferentially in the gas-rich nuclear regions surrounding active galactic nuclei (AGNs; McKernan et al. 2018)—as a result, e.g., of Kozai–Lidov resonances (Antonini & Perets 2012) or simply fragmentation/instabilities of the AGN accretion disk (Stone et al. 2017). Furthermore, accretion onto stellar-origin or intermediate-mass black holes has been proposed as an explanation for ultraluminous X-ray sources (see, e.g., Miller et al. 2004). Accretion, in combination with mergers, is also thought to be the main channel via which black-hole seeds evolve into the supermassive black holes we observe today.

Therefore, at least some SOBHBs or IMBHBs may still be accreting matter in the LISA band and perhaps even at merger. The accretion-driven EM emission may not have been detected because these sources are too far,<fnref n="10" /> because accretion is radiatively inefficient (Frank et al. 2002), or because the sky position uncertainty provided by GWs is too large for follow-up campaigns. Note also that LISA is expected to detect up to several tens of SOBHBs (Sesana 2016; Tamanini et al. 2020). If only one such system were accreting, and if the possibility for accretion were not included in the GW templates used for the analysis, the parameter estimation may mistakenly point toward a modification of general relativity (Barausse et al. 2016; Carson & Yagi 2020a; Gnocchi et al. 2019)—a claim that would have groundbreaking effects on physics. Furthermore, LISA may provide an accurate sky localization for these sources, thus increasing the chances of detecting a putative EM counterpart, with important implications for multi-messenger astronomy and cosmology.

With these motivations, in this work, we analyze the effect of gas accretion on standalone IMBHB LISA detections and on joint LISA+ground multiwavelength SOBHB observations. We find that accretion introduces a −4 Post-Newtonian (PN) correction to the phase,<fnref n="11" /> thus potentially dominating over the GW-driven evolution at low frequencies (see also, e.g., Holgado & Ricker 2019). The systems we consider will be driven by gravitational-wave emission, with accretion acting as a perturbative correction and, therefore, leaving an imprint on the GW phasing. We explore the consequences of this fact for GW parameter estimation, i.e., we assess both with what uncertainty the accretion rate can be recovered when the possibility for accretion is included in the templates and how much the estimate of the binary parameters will be biased if it is not. We also look at the prospects of identifying the EM

emission from accreting SOBHBs and IMBHBs with observational facilities available when LISA flies.

In Section 2, we begin by summarizing the effect of accretion on the GW waveform and on the binary evolution. In Section 3, we describe how we generate astrophysical catalogs and simulate future detections. We present our results in Section 4, and we summarize them in Section 5. We use geometrized units in which $G = c = 1$. We denote the total mass by $M = m_1 + m_2$, the reduced mass by $\mu = m_1m_2/M$, and the chirp mass by $M = (\mu^3 M^2)^{1/5}$.

## 2. Shift of the Merger Time and Waveform Corrections Due to Accretion

Let us parameterize the mass accretion rate of each component of a (circular) black-hole binary (with masses $m_i$, $i = 1, 2$) by the Eddington ratio

$$f_{\text{Edd},i} = \frac{\dot{m}_i}{\dot{m}_{\text{Edd}}}, \qquad (1)$$

where $\dot{m}_{\text{Edd}} \approx 2.2 \times 10^{-8} \left(\frac{m_i}{M_{\odot}}\right) M_{\odot} \text{yr}^{-1}$ is the Eddington accretion rate (obtained from the Eddington luminosity assuming radiative efficiency $\eta = 0.1$). Since the accretion timescale exceeds the dynamical timescales of the binary when the latter is in the frequency band of LISA or ground detectors, the effect on the phase can be computed using the stationary-phase approximation, and to the leading order at low frequencies it reads (see derivation in Appendix A)

$$\phi_{\text{acc}} \sim \alpha f_{\text{Edd}} f^{-13/3}, \qquad (2)$$

where $f$ is the GW frequency and $\alpha$ is a coefficient that depends on the binary parameters. Since the leading-order term in the phase in vacuum is $\sim f^{-5/3}$ (Maggiore 2008), this is a −4PN term, which dominates the binary evolution at low frequencies. In the frequency range of LISA observations, due to the smallness of the prefactor, this term will be a small correction to the vacuum GW phase. In other words, our SOBHB and IMBHB sources will emit GWs well above the frequency at which accretion becomes subdominant,

$$f \gg 1.1 \times 10^{-4} \left(\frac{f_{\text{Edd}}}{1}\right)^{3/8} \left(\frac{M}{10 M_{\odot}}\right)^{-5/8} \text{Hz}, \qquad (3)$$

see Appendix A.

As a result of accretion, the phase evolution accelerates and the binary merges earlier (i.e., in less time and in fewer GW cycles) than in vacuum. Note that in this work, we neglect for simplicity the hydrodynamic drag produced by the transfer of linear momentum by the accreting gas (Barausse & Rezzolla 2008; Chen & Shen 2019), which would further contribute to the shift of the merger time (see Appendix A).

In Figure 1, we show the time $T$ needed for an SOBHB to enter the band of ground detectors (top panel), the time difference $\Delta T$ in the merger time induced by accretion (middle panel), and the difference $\Delta\phi$ in the total (accumulated) GW phase due to accretion (bottom panel), as functions of the initial GW frequency in the LISA band and for various SOBHB masses. All of these quantities can be computed either numerically solving Equation (18) or using the perturbative expansions in Equations (22)–(23). The two approaches are in

<footnote n="10">¹⁰ Note for instance that accreting black holes in X-ray binaries are mostly observed in the Galaxy, with only a few observed in nearby galaxies. Among the latter, the farthest is M83 (Ducci et al. 2013), which is only ~4.5 Mpc away, versus the several hundred Mpc of the LIGO/Virgo SOBHBs (Abbott et al. 2019).</footnote>

<footnote n="11">¹¹ In the GW phase, an $n$PN correction scales as $(v/c)^{2n} \sim f^{2n/3}$ (v and $f$ being the binary’s orbital velocity and the GW frequency) relative to the leading-order general relativistic term.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Time needed to enter the band of terrestrial detectors neglecting accretion (top panel), time difference caused by accretion (middle panel), and corresponding GW phase difference (bottom panel), as functions of the initial GW frequency for three equal-mass SOBHBs. We choose $f_{Edd} = 1$ as a reference, since the time and phase differences scale linearly with $f_{Edd}$. The full (empty) circles mark the points corresponding to $T = 10$ yr ($T = 4$ yr). The systems to the right of the full (empty) circles therefore have $T < 10$ yr ($T < 4$ yr).</figcaption></figure>

excellent agreement because the contribution of accretion is subdominant in all cases.

As a useful rule of thumb, time differences $\Delta T > 10$ s (Sesana 2016) and phase differences $\gtrsim 1$ rad (Flanagan & Hughes 1998; Lindblom et al. 2008) are large enough to be detectable.

For low initial frequency, the effect of accretion on $\Delta T$ and on the phase is stronger, but the time $T$ is also very large, i.e., multiwavelength observations will be impossible in practice. One may try to detect accretion with LISA data alone, but note that the mission's duration will not exceed 10 yr (with a nominal duration of 4 yr), due to the finite consumables carried by the spacecraft. For these reasons, we mark in Figure 1 the phase and time differences for an SOBHB that enters the band of ground detectors in 10 (4) yr by full (empty) circles. The part of the

curves to the right of these circles then corresponds to $T < 10$ yr ($T < 4$ yr), which would make a joint LISA+ground detection possible in practical terms. Overall, the results of Figure 1 (which scale linearly with $f_{Edd}$) suggest that only $f_{Edd} > 0.1$ would give a potentially detectable effect, i.e., $\Delta T > 10$ s and $\Delta\phi \gtrsim 1$. We will verify this with more rigorous techniques in the following.

### 3. Measuring Accretion Effects for SOBHBs and IMBHBs

In order to quantify the ability of multiband SOBHB detections and standalone IMBHB observations to constrain the accretion model, we perform two analyses: (i) a simple Fisher-matrix analysis to explore the whole parameter space, and (ii) a more refined Monte Carlo Markov Chain (MCMC) analysis for the best-candidate events. Note that the Fisher-matrix analysis is only valid for large signal-to-noise ratios (S/Ns) (Vallisneri 2008). Therefore, we expect it to provide only qualitatively correct results for SOBHBs in the LISA band (for which the S/N is at most 15–20 in the most optimistic cases, see below). Nevertheless, we expect the Fisher-matrix analysis to be accurate for the IMBHBs we consider, for which S/N = $\mathcal{O}(100)$.

In both the Fisher and MCMC analyses, we only account for the contribution due to accretion in the GW phase, and neglect the subleading contribution to the amplitude. Since accretion is important at low frequency, high-order PN terms (including the spin) should be irrelevant for our analysis, but we include them for completeness and to estimate possible correlations.

For simplicity, in the Fisher analysis, we also neglect the motion of the antenna during the observation. This is instead included in the MCMC analysis, in order to estimate the ability to localize the source in the sky and measure the accretion rate at the same time.

Finally, we consider two situations: one (referred to as LISA+Earth) in which we simulate a multiband SOBHB detection (LISA combined with a ground-based interferometer) and another (referred to as LISA-only) in which we simulate a standalone (either SOBHB or IMBHB) detection by LISA. In the LISA+Earth case, to simulate a multiband detection, one can follow two options: combine statistically the noise curves of LISA with that of a given ground-based detector or, alternatively (but less rigorously), assume that the merger time can be computed independently by the ground-based detector, so that the dimension of the parameter space of the analysis is effectively reduced. In the Fisher analysis, we follow the latter, simpler approach, and we therefore effectively remove the merger time from the template parameters in the LISA+Earth case. In the MCMC analysis, we keep $t_c$ as a free parameter, restricting it by using a narrow prior. In all cases, we adopt the LISA noise curve reported by Audley et al. (2017), whose high-frequency part is based on a single-link optical measurement system noise of $10 \text{ pm}/\sqrt{\text{Hz}}$.

#### 3.1. Fisher Analysis and Event Rates for SOBHBs

In the Fisher analysis, we adopt a TaylorF2 template approximate for spinning binaries up to 3.5PN order (Droz et al. 1999), with the addition of the leading-order accretion term presented in Equation (2). Therefore, our GW template for the Fisher analysis has seven parameters (masses, merger time and phase, the two dimensionless spins $\chi_{1,2}$, besides the Eddington accretion ratio $f_{Edd}$).

Given a waveform template $h(\zeta, f)$ in the frequency domain and a set of waveform parameters $\zeta$, the error associated with the
measurement of parameter $\zeta^a$ (with all other parameters margin-
alized upon) is $\sigma_a = \sqrt{\sum^{aa}}$, where the covariance matrix $\Sigma^{ab}$ is
given by the inverse of the Fisher matrix, $\Gamma_{ab} = (\partial_{\zeta^a}h|\partial_{\zeta^b}h)_{\zeta=\zeta_0}$.
Here, $\zeta_0$ are the injected values of the parameters, and the inner
product is defined by

$$ (g|h) = 4 \operatorname{Re} \int_{f_0}^{f_{\max}} df \frac{\tilde{h}(f) \tilde{g}^*(f)}{S_h(f)}, \quad (4) $$

where $S_h(f)$ is the detector noise spectral density.

While the number (and the very existence) of IMBHBs in the LISA band is very uncertain, our Fisher-matrix analysis, coupled with simulated astrophysical populations calibrated to the LIGO/Virgo data, can easily provide estimates of the number of SOBHBs detectable by LISA for which accretion can be measured. The intrinsic number of SOBHBs merging per (detector-frame) unit time and (source-frame) masses is given by (Hartwig et al. 2016)

$$ \frac{d\dot{N}}{dm_1 dm_2} = \int dz R \frac{dt}{dz} \frac{d^2 p}{dm_1 dm_2} 4\pi d_C^2, \quad (5) $$

where $d_C$ is the comoving distance, $R = 53.2 \text{ Gpc}^{-3} \text{ yr}^{-1}$ is the best estimate for the intrinsic merger rate measured by the first and second LIGO/Virgo runs (Abbott et al. 2019), the probability distribution function for the source-frame masses —$d^2p/dm_1dm_2$— is given by “model B” of Abbott et al. (2019), while

$$ \frac{dt}{dz} = \frac{1}{H_0(1+z)\sqrt{\Omega_m(1+z)^3 + \Omega_\Lambda}} $$

is computed using our fiducial cosmology $H_0 = 67.9 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.306$, $\Omega_\Lambda = 0.694$ (Ade et al. 2016). In order to obtain synthetic astrophysical catalogs of merging as well as inspiraling sources, we use Equation (5) to simulate mergers in a period much longer than the LISA mission duration, by assuming a uniformly distributed merger time $t_c$. The latter can be easily converted into the initial GW frequency $f_0 = [5/(256 t_c)]^{3/8} M^{-5/8}/\pi$, where $f_0$, $t_c$, and the chirp mass $M$ must be computed in the same (detector- or source-) frame.

We constrain the comoving distance in the range $d_C \le 2 \text{ Gpc}$ and the initial source-frame GW frequency in the range $f_0 \in [4 \text{ mHz}, 10 \text{ Hz}]$. For the chosen mass model, we generate 20 realizations, and for each realization, we consider two LISA mission durations (4 or 10 yr), for a total of 40 catalogs.

In the LISA-only case for SOBHBs, we assume that a single event within the catalog is detected if either of the following conditions occurs (Moore et al. 2019; Tamanini et al. 2020):

$t_c < 100 \text{ yr and } S/N \ge 15$, or

$t_c > 100 \text{ yr and } S/N \ge 10,$

where the latter S/N threshold is lower because binaries with long merger times are accurately described by Newtonian waveforms in the LISA band (Mangiagli et al. 2019) and can be therefore detected by a different search strategy (Tamanini et al. 2020), akin, e.g., to the one used for white-dwarf binaries.

In the LISA+Earth case, the S/N threshold is lower for events that can be detected on Earth (Moore et al. 2019; Tamanini et al. 2020)

$$ t_c < 10 \text{ yr and } S/N \ge 9.5. \quad (6) $$

These events would indeed be detected through an archival search following their ground-based detection.

## 3.2. MCMC with Sky Localization and Antenna Motion

For the MCMC analysis, we adopt the PhenomD template (Husa et al. 2016; Khan et al. 2016) with the inclusion of the phase term due to accretion. In this case, we also account for the motion of the antenna during the observation. More specifically, the standard part of the GW template is the same as in Tamanini et al. (2020) and contains five additional parameters besides those adopted for the Fisher-matrix analysis: two angles identifying the source position with respect to the detector ($\phi, \theta$), the GW polarization ($\psi$), the inclination of the system ($i$), and the luminosity distance ($d_L$).

In the LISA-only scenario, we use $f_0$ as sampling parameter and assume a flat prior for it. In the LISA+Earth scenario, we use $t_c$ with a Gaussian prior centered around the true value with width $\sigma_{t_c} = 10^{-3}$ s, which models the fact that $t_c$ can be measured with great precision in this scenario. For IMBHBs, we consider a single LISA-only scenario.

When including the source location, different realizations of the angles for the same astrophysical system yield different S/Ns. This affects the precision within which one can recover the parameters of the source, including the sky position itself and $f_{\text{Edd}}$. In order to cross-check results obtained with our Fisher-matrix analysis and to quantify this variability, we select from the catalog an astrophysical system for which the accretion parameter can be measured precisely through the Fisher-matrix approach, and draw three different realizations of ($\phi, \theta, \psi, i$) yielding a low S/N ~ 9, a medium S/N ~ 15, and a high S/N ~ 20, respectively. The medium S/N system is chosen so that its S/N is close to the value obtained by averaging over the angles.

For IMBHBs, we consider two systems (see details in the next section): one merging in the LIGO/Virgo band, and one with higher masses, merging at lower frequencies. We choose the initial frequency so that both systems merge in 10 yr, the longest possible LISA mission duration.

For each of these systems, we perform a full Bayesian analysis (see Appendix B). We simulate GW data $d(f)$ as it would be measured by LISA, computing the response of the detector (accounting for the constellation's motion) by following Marsat & Baker (2018). We work in the zero-noise approximation in order to speed up the computation. Adding noise to the GW signal should not affect the parameter estimation drastically, leading mostly to a displacement in the maximum of the parameter distribution (Rodriguez et al. 2014).

We perform two different analyses: in the first one, we generate data with a nonzero value for $f_{\text{Edd}}$ and include it as a free parameter in the Bayesian analysis, in order to estimate with what precision it can be recovered. In the second case, data are also generated with a nonzero value of $f_{\text{Edd}}$, but when doing the analysis, we set $f_{\text{Edd}} = 0$ in the templates, in order to measure the bias in the parameter estimation. In all cases, the posterior distribution is computed using Bayes' theorem. Additional details are given in Appendix B.

# 4. Results

## 4.1. Event Rates for SOBHBs

For the simulated astrophysical populations, we first use a Fisher-matrix analysis to quantify the possibility to measure
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Marginalized distributions of $f_{Edd}$ for the considered SOBHB for various realizations of the angles, and in the scenarios LISA+Earth (green) and LISA-only (red). Red lines indicate the injected value of the accretion rate, $f_{Edd} = 1$. In the LISA+Earth scenario and for higher S/N the marginalized distribution is strongly peaked but still consistent with $f_{Edd} = 0$.</figcaption></figure>

<table><caption>**Table 1**
Number of Detectable SOBHB Events for Various Configurations</caption><thead><tr><th rowspan="2">Duration</th><th colspan="5">LISA+Earth</th><th colspan="5">LISA-only</th></tr><tr><th>All</th><th>f<sub>Edd</sub></th><th>100%</th><th>50%</th><th>10%</th><th>All</th><th>f<sub>Edd</sub></th><th>100%</th><th>50%</th><th>10%</th></tr></thead><tbody><tr><td rowspan="2">4 yr</td><td rowspan="2">88 ± 8</td><td>1</td><td>0.1 ± 0.2</td><td>0</td><td>0</td><td rowspan="2">77 ± 8</td><td>1</td><td>0</td><td>0</td><td>0</td></tr><tr><td>10</td><td>4.1 ± 2.3</td><td>1.7 ± 1.2</td><td>0.1 ± 0.2</td><td>10</td><td>1.6 ± 1.4</td><td>0.6 ± 0.6</td><td>0</td></tr><tr><td rowspan="2">10 yr</td><td rowspan="2">207 ± 11</td><td>1</td><td>5.2 ± 1.9</td><td>1.1 ± 1.2</td><td>0.1 ± 0.2</td><td rowspan="2">182 ± 10</td><td>1</td><td>1.5 ± 1.2</td><td>0.4 ± 0.7</td><td>0</td></tr><tr><td>10</td><td>36 ± 4</td><td>32 ± 3</td><td>5.2 ± 1.9</td><td>10</td><td>11 ± 3</td><td>9.5 ± 2.7</td><td>1.5 ± 1.2</td></tr></tbody></table>

Note. “All” stands for the total number of detectable events, whereas 100%, 50%, and 10% stand for the number of events for which $f_{Edd}$ is measured with a relative error of 100%, 50%, and 10%, respectively, according to the Fisher analysis. All numbers are averaged over 20 catalogs and presented with $1\sigma$ errors. Super-Eddington accretion will be detectable for a good fraction of multiband events if the LISA mission duration is 10 yr.

$f_{Edd}$ at a given precision. Table 1 shows the average number of detected SOBHBs, and the number of SOBHBs for which $f_{Edd}$ can be measured within a given precision. The results are obtained by averaging the Fisher matrix over sky position and source inclination (while neglecting, as already mentioned, the LISA constellation’s motion), for different injected values of $f_{Edd}$. Our results for the total number of detected events are consistent with Tamanini et al. (2020) and Sesana (2016).

In particular, for the LISA+Earth case and a 10 yr mission, super-Eddington accretion $f_{Edd} \approx 10$ can be measured within 50% precision in about 15% of the total detectable events ($\approx 200$), while a measurement within 10% is only possible in ~2% of the events. Note that the statistical errors scale approximately linearly with $f_{Edd}$. Therefore, when injecting a lower accretion rate, the number of events for which accretion is measurable is significantly smaller. For example, $f_{Edd} = 1$ is marginally detectable in $\lesssim 1$ event in the most optimistic scenario, whereas smaller values of the accretion rates are not measurable.

As expected, a multiband observation improves the measurements of a negative-PN term, including the -4PN term due to accretion: the event rates for the LISA-only case are thus smaller by a factor of a few relative to the LISA+Earth case.

## 4.2. Measuring Accretion and Sky Localization

For our MCMC analysis, we select one representative SOBHB system from our synthetic astrophysical catalogs, and choose two optimistic IMBHB systems on the basis of a Fisher-matrix analysis spanning the parameter space, i.e., the errors on $f_{Edd}$ provided by the chosen IMBHBs are roughly the smallest throughout the parameter space. In more detail, the systems that we consider are

1. An SOBHB with $m_1 = 42.1 M_\odot$, $m_2 = 39.8 M_\odot$, $\chi_1 = 0.008$, $\chi_2 = 0.44$, at a distance $d_L = 416$ Mpc;

2. An IMBHB with $m_1 = 315 M_\odot$, $m_2 = 284 M_\odot$, $\chi_1 = 0.9$, $\chi_2 = 0.85$, referred to as “light IMBHB;”

3. Another IMBHB with $m_1 = 1000 M_\odot$, $m_2 = 900 M_\odot$, $\chi_1 = 0.9$, $\chi_2 = 0.85$, referred to as “heavy IMBHB.”

For all three sources, we set $t_c \approx 10$ yr. We study the IMBHB systems at two different redshifts, $z = 0.1$ and $z = 0.5$, in order to estimate up to what distance the presence of accretion in the binary would be detectable. The IMBHBs’ masses are in the source frame and are kept fixed when the redshift is changed.

For each realization of the angles $(\theta, \phi, \iota, \psi)$, we compute and sample the posterior distribution as explained in Appendix B. As expected, the precision of the parameter measurements increases
<figure><img src="image_4.png" /><figcaption>Figure 3. Same as Figure 2 but for an injected value of $f_{Edd} = 10$. Accretion is detected in all scenarios and for any angle realization.</figcaption></figure>

with the S/N. We find that the accretion parameter is strongly correlated with the intrinsic parameters of the source ($M, \mu/M, f_0, \chi_s, \chi_a$), where $\chi_s$ and $\chi_a$ are defined in Appendix B.

In Figure 2, we show the marginalized distributions of $f_{Edd}$ for the chosen SOBHB system, for various S/Ns and for an injected value of $f_{Edd} = 1$.

For the high- and the medium-S/N cases, already in the LISA-only scenario, the posteriors indicate the presence of accretion. The marginalized distribution for $f_{Edd}$ can be compared with those obtained when constraining modifications of GR (some of which affect the vacuum waveform in a similar fashion as accretion, i.e., at negative PN orders) in the parameterized post Einsteinian framework (Yunes & Pretorius 2009). In that case, as discussed in an upcoming paper (A. Toubiana et al. 2020b, in preparation), the marginalized distribution of the non GR-parameters is mostly flat up to a threshold (representing the upper bound that can be placed on the parameters under scrutiny), and then goes to 0. In contrast, we see in Figure 2 that for high and medium S/N in the LISA-only scenario, the distribution peaks at some nonvanishing value, indicating the presence of a nonzero modification to the vacuum waveform.

In Figure 3, we show the same as in Figure 2, but for an injected value of $f_{Edd} = 10$. This high accretion rate can be detected more easily even in the low-S/N case and in the LISA-only scenario, since in this case, $f_{Edd} = 0$ is outside the support of the distribution. Thus, for super-Eddington accreting binaries in the LISA band, there is a concrete chance to detect the effect of high rates of accretion on the waveform for most SOBHB events.

In Table 2, we show the recovered 68% confidence intervals (CI) and median values for $f_{Edd}$ and the sky localization ($\Delta\Omega$). In the $f_{Edd} = 1$ case, since the distribution is leaning against the boundary of the prior (see Figure 2), we define the 68% CI for $f_{Edd}$ by taking the lower 68% values. Instead, in the $f_{Edd} = 10$ case, the interval is centered around the median values. The marginalized distributions for $\Delta\Omega$ are approximately Gaussian and are centered around the injected value. Thus, we define the solid angle as (Cutler 1998):

$$ \Delta\Omega = 2\pi\sqrt{(\Sigma^{\phi,\phi})(\Sigma^{\cos(\theta),\cos(\theta)} - (\Sigma^{\phi,\cos(\theta)})^2)}. \quad (7) $$

<caption>Table 2</caption>

Recovered 68% CI on the Accretion Parameter $f_{Edd}$ and on the Sky Localization $\Delta\Omega$ for SOBHBs and for Various Realizations in the LISA+Earth Scenario

<table><thead><tr><th rowspan="2"></th><th colspan="2">$f_{Edd}^{\text{injected}} = 1$</th><th colspan="2">$f_{Edd}^{\text{injected}} = 10$</th></tr><tr><th>$f_{Edd}$</th><th>$\Delta\Omega$ (deg<sup>2</sup>)</th><th>$f_{Edd}$</th><th>$\Delta\Omega$ (deg<sup>2</sup>)</th></tr></thead><tbody><tr><td>High S/N</td><td>0.68<sup>+0.02</sup><sub>-0.68</sub></td><td>0.14</td><td>9.46<sup>+0.53</sup><sub>-0.83</sub></td><td>0.14</td></tr><tr><td>Medium S/N</td><td>0.70<sup>+0.25</sup><sub>-0.70</sub></td><td>0.06</td><td>9.23<sup>+0.77</sup><sub>-1.23</sub></td><td>0.06</td></tr><tr><td>Low S/N</td><td>1.18<sup>+0.50</sup><sub>-1.18</sub></td><td>0.33</td><td>8.82<sup>+1.85</sup><sub>-2.96</sub></td><td>0.34</td></tr><tr><td>Fisher matrix</td><td>1.00<sup>+1.20</sup><sub>-1.20</sub></td><td>...</td><td>10.00<sup>+1.20</sup><sub>-1.20</sub></td><td>...</td></tr></tbody></table>

Note. The last row gives the statistical error estimated with a Fisher-matrix analysis. The presence of accretion should be detected for super-Eddington accreting systems. The error on the sky position is always within the field of view of Athena and SKA, allowing (potentially) for electromagnetic follow-up.

<caption>Table 3</caption>

Recovered 68% CI on the Accretion Parameter $f_{Edd}$ and on the Sky Localization $\Delta\Omega$ for the IMBHBs Considered in This Work, at Redshift $z = 0.1$

<table><thead><tr><th rowspan="2"></th><th colspan="2">$f_{Edd}^{\text{injected}} = 0.1$</th><th colspan="2">$f_{Edd}^{\text{injected}} = 1$</th></tr><tr><th>$f_{Edd}$</th><th>$\Delta\Omega$ (deg<sup>2</sup>)</th><th>$f_{Edd}$</th><th>$\Delta\Omega$ (deg<sup>2</sup>)</th></tr></thead><tbody><tr><td>Light IMBHB</td><td>0.14<sup>+0.04</sup><sub>-0.14</sub></td><td>0.01</td><td>1.05<sup>+0.10</sup><sub>-0.11</sub></td><td>0.01</td></tr><tr><td>Heavy IMBHB</td><td>0.17<sup>+0.05</sup><sub>-0.17</sub></td><td>0.007</td><td>1.04<sup>+0.15</sup><sub>-0.27</sub></td><td>0.006</td></tr><tr><td>Fisher matrix</td><td>0.10<sup>+0.38</sup><sub>-0.38</sub></td><td>...</td><td>1.00<sup>+0.38</sup><sub>-0.38</sub></td><td>...</td></tr></tbody></table>

Note. The statistical error estimated with our Fisher-matrix analysis is similar for the two IMBHBs. The presence of accretion should be detected for Eddington accreting systems. As for SOBHBs, the error on the sky position is always within the field of view of Athena and SKA.

We show the same quantities for our IMBHB events in Table 3. There, in the case $f_{Edd} = 1$, we define the 68% CI for $f_{Edd}$ centered on the median, and in the case $f_{Edd} = 0.1$, we define it by taking the lower 68% values. In all cases considered here, the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Marginalized distributions of $f_{Edd}$ for our two IMBHB systems at redshifts $z = 0.1$ (red) and $z = 0.5$ (green) for an injected value of $f_{Edd} = 1$. Accretion can be measured in both systems at $z = 0.1$, but not at higher redshift.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Marginalized distributions of $f_{Edd}$ for our two IMBHB sources for injected values of $f_{Edd} = 0.1$ (green) and $f_{Edd} = 1$ (red) at $z = 0.1$. Accretion at this redshift needs to be approximately Eddington-level or stronger to be measured.</figcaption></figure>

error on the sky localization is much smaller than the nominal field of view of future X-ray and radio missions, potentially allowing for the detection of electromagnetic counterparts. We will discuss this possibility in Section 4.3.

While overall in qualitative agreement, the differences between Fisher-matrix and MCMC results could be due to the effect of the priors, to the non-Gaussianity of the posterior distribution, to the treatment of the angles, and/or to the finite S/N of the sources considered. Nonetheless, the predicted errors on $f_{Edd}$ are of the same order of magnitude in both treatments, confirming the main conclusions we drew for SOBHBs using the Fisher analysis.

In Figure 4, we compare how well can we recover $f_{Edd}$ for IMBHBs at different redshifts, for injected $f_{Edd} = 1$. If the system is too far, the distribution tends to be flat and the effect of accretion is hardly noticeable. This is because of the lower S/N, but also because the detector-frame mass becomes larger

at higher redshift, speeding up the evolution of the system and thus providing less information on negative PN-order modifications.

Finally, in Figure 5, we show how well can we recover $f_{Edd}$ in IMBHBs for an injected values of $f_{Edd} = 0.1$ at $z = 0.1$. As in the case of SOBHBs commented above, the marginalized distribution is compatible with $f_{Edd} = 0$, but the presence of a clear peak at $f_{Edd} \neq 0$ favors the presence of accretion.

### 4.2.1. Estimating Biases

The above results indicate that if accretion is present, it could lead to a measurable change in the GW signal. Thus, if accretion is not taken into account, the estimation of other source parameters could be significantly biased. Since $f_{Edd}$ correlates mostly with the intrinsic parameters of the source, the latter should be the most affected.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Bias in the SOBH binary masses and time to coalescence induced by ignoring the corrections due to accretion when $f_{Edd} = 1$, for various angle realizations in the *LISA-only* scenario (blue) compared to the displacement found in vacuum systems. Boxes and whiskers delimit the 50% CI and the 90% CI, respectively, and both are centered around the median, indicated by lines inside the boxes. For this level of accretion, bias is not significant, even for the high-S/N realization.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same as in Figure 6 but comparing $f_{Edd} = 0$ with $f_{Edd} = 10$. Bias here is significant for all S/N realizations.</figcaption></figure>

For SOBHBs in the *LISA-only* scenario, we find that in the three cases (high, medium, and low S/N), the signal can be recovered by an effectual template with $f_{Edd} = 0$, i.e., we find a maximum for the posterior distribution that, in the worst cases, can be incompatible with the injected real value. The S/N of this effectual template is very similar to the injection's S/N ($S/N_{inj} - S/N_{eff} \lesssim 0.7$), and could thus trigger a detection. The bias in the parameter estimation and the relative drop in S/N is higher for lower-S/N systems and for higher injected accretion rates. The effectual template, in particular, has a higher chirp mass and a higher mass ratio, while the initial frequency is shifted toward higher values. In Figures 6 and 7, we show how this impacts the estimate of the masses and time to coalescence for two representative values, $f_{Edd} = 1$ and $f_{Edd} = 10$. In both cases, we compare to the recovered distribution of masses for vacuum GR. The mass of the primary black hole is shifted toward higher values, whereas the secondary mass gets lower. As a result, the time to coalescence is underestimated. For super-Eddington accretion, this shift in time to coalescence is at the level of tens of seconds. A multiband observation could then help identify a bias

due to accretion in the parameter estimation, since ground-based detectors would measure very precisely the time to coalescence when the signal enters in their band (Sesana 2016).

In order to estimate this possibility, we repeat the above analysis in the *LISA+Earth* scenario. In this case, the time to coalescence is constrained to within 1 ms from its true value, so no bias in $t_c$ is possible. Nevertheless, signals can still be recovered by an effectual template, although, with a larger mismatch from the true signal. In Figure 8, we show the difference between the recovered masses and total mass and the injected values. More in general, ground-based detectors should be largely insensitive to these low-frequency terms, as discussed in Carson & Yagi (2020b). In this forecast study, the projected constraint on -4PN terms with the planned third-generation detector Cosmic Explorer is ten orders of magnitude worse than the projected constraint with LISA. We thus expect that observations with ground-based detectors should not be biased by omitting -4 PN terms. Therefore, for values of $f_{Edd}$ for which the LISA parameter estimation is significantly biased, the posterior distributions obtained with LISA and with
<figure><img src="image_9.png" /><figcaption>Figure 8. Same as in Figure 6 but for the LISA+Earth scenario. We show $m_{tot}$ rather of $t_c$, the latter being fixed by the narrow prior in this scenario. Bias in the masses can still be significant for medium and high S/N realizations, despite the constraint provided by ground-based detectors.</figcaption></figure>

ground-based detectors might not even be compatible, which would hint at an unmodeled effect.

It is noteworthy that for the SOBHB events, the sky localization is barely affected by accretion and remains excellent, as the distribution remains a Gaussian centered around the injected values with errors similar to the ones shown in Table 2. In the case of IMBHBs, on the other hand, there is also a bias in the sky localization, i.e., the injected value may lie outside the 90% CI. This is due to the very small errors in sky position, and in fact, the true localization is very close to the recovered one, within 0.05 deg². Therefore, for most realistic purposes, the sky localization is satisfactorily recovered.

Since we did not consider any modification to the GW amplitude, there is no strong correlation between $f_{Edd}$ and the luminosity distance $d_L$. Thus, when fixing $f_{Edd} = 0$ as we did here, there is no bias on the estimation of $d_L$, contrary to the Fisher-matrix analysis in Tamanini et al. (2020), who also used waveforms modifying GR at -4PN order in phase but included the leading-order modification to the amplitude too.

## 4.3. Prospects for Multiband and Multi-messenger Astronomy

According to our MCMC analysis, both SOBHBs and IMBHBs can be localized in the sky to within the fields of view of X-ray and radio instruments such as Athena WFI and SKA, $\Delta\Omega_{Athena} = 0.4$ deg², $\Delta\Omega_{SKA} = 0.5$ deg² (Dewdney et al. 2014; Meidinger 2018). This will allow the relevant region of the sky to be covered in a single viewing,<fnref n="12" /> thus, potentially allowing for the coincident detection of an X-ray and/or radio counterpart to strongly accreting black-hole binaries. Even if the sky localization was biased, as might be the case for IMBHBs, we estimated that the true position would still fall inside the field of view of the instruments. In the following, we compute the X-ray and radio emission of the binaries, and estimate the necessary integration time for detection by a single instrument viewing.

We start by estimating the X-ray flux. For this purpose, we assume that the accretion process has radiative efficiency $\eta = 0.1$ (which is a good approximation at $f_{Edd} < 1$) and that only a fraction $\eta_X = 0.1$ of the EM radiation is emitted in X-rays (“bolometric correction”). We find the X-ray flux from a single accreting black hole to be

$$F_X \approx 1 \times 10^{-13} f_{Edd} \left(\frac{M}{M_{\odot}}\right) \left(\frac{\text{Mpc}}{d_L}\right)^2 \text{erg cm}^{-2} \text{s}^{-1}. \quad (8)$$

This should be compared with the flux sensitivity of the Athena WFI for a given integration time, $T_{int}$. Following McGee et al. (2020), Athena’s flux sensitivity for a 5σ detection is

$$F_X^{\text{Athena}} = 1 \times 10^{-15} \left(\frac{10^3 \text{s}}{T_{\text{int}}}\right)^{1/2} \text{erg cm}^{-2} \text{s}^{-1}. \quad (9)$$

The minimum integration time for a binary where only one black hole is emitting is then given by

$$T_{\text{int}} \approx 8 \times 10^{-2} f_{\text{Edd}}^{-2} \left(\frac{d_L}{\text{Mpc}}\right)^4 \left(\frac{M_{\odot}}{M}\right)^2 \text{s}. \quad (10)$$

Note that if the two black holes have similar mass and are both accreting, the cumulative flux is given by twice the value in Equation (8), and therefore, the minimum integration time is one-fourth of that in Equation (10).

For the best-candidate SOBHB event in our synthetic astrophysical catalogs, the required exposure time is $T_{int} \gtrsim 1 \times 10^6 f_{Edd}^{-2}$ s. Thus, even if we were to assume $f_{Edd} \approx 1$, the integration time would have to be of several days. Assuming super-Eddington accretion $f_{Edd} > 1$ is unlikely to help as the radiative efficiency is expected to be considerably lower than our assumed $\eta = 0.1$, i.e., the bolometric luminosity is not expected to significantly exceed the Eddington luminosity (Shakura & Sunyaev 1973; Poutanen et al. 2007; Sadowski 2011). Moreover, as previously discussed, high accretion rates in SOBHBs likely require environments with large gas densities, whose optical thickness further reduces the chances of an EM detection. For the considered IMBHB systems, the required integration time is between 24 and 2 hr for Eddington-level accretion, for the light and heavy systems, respectively. This estimate suggests that

<footnote n="12">¹² In some cases, the correlation between the sky position angles can imprint an asymmetric shape to the localized region, which might therefore partially fall outside of the field of view. However, this would still only require $\mathcal{O}(1)$ viewings.</footnote>
detection of X-ray counterparts will be possible for highly accreting IMBHs.

A binary system in external magnetic fields may also launch dual radio jets, which get amplified by the coalescence (Palenzuela et al. 2010) relative to similar jets observed in isolated black holes (Steiner et al. 2012). See also Moesta et al. (2012) for simulations that yield ~100 times larger (though less collimated) fluxes than Palenzuela et al. (2010). Assuming a fiducial value $\eta = 0.1$ for the radiative efficiency of the process and $\eta_{\text{radio}} = 0.1$ for the fraction of emission in the radio band, the corresponding peak flux<fnref n="13" /> is (Palenzuela et al. 2010; Tamanini et al. 2016)

$$F_{\text{flare}} \approx 2 \times 10^{-13} f_{\text{Edd}} q^2 \left(\frac{D_L}{\text{Mpc}}\right)^{-2} \left(\frac{M}{M_{\odot}}\right) \times \text{erg cm}^{-2} \text{s}^{-1}, \quad (11)$$

where $q \le 1$ is the mass ratio. The flare flux can then be compared with the SKA-mid sensitivity in the phase 1 implementation. The required sensitivity at frequency $\nu_{\text{SKA}}$ for SKA,

$$F_{\text{SKA}} = 5 \times 10^{-16} \left(\frac{10^{-2} \text{s}}{T_{\text{obs}}}\right)^{1/2} \left(\frac{\nu_{\text{SKA}}}{\text{GHz}}\right) \text{erg cm}^{-2} \text{s}^{-1}, \quad (12)$$

is reached for an observation time $T_{\text{obs}} \sim 10^{-2} \text{s}$ for our best SOBHB event. The observation time should be smaller than the duration of the merger (i.e., the duration of the flare) for the system (Steiner et al. 2012), $T_{\text{flare}} \sim 25 \frac{M}{100 M_{\odot}} \text{ms}$. This condition is not satisfied for SOBHBs. There is, however, the concrete possibility to detect a signal in the radio band for IMBHs, for which for the light and heavy systems $T_{\text{obs}} \approx 40 - 4 \text{ms} < T_{\text{flare}}$. The performance of full SKA should improve by an order of magnitude with respect to Equation (12), reducing the required integration time by a factor 100.

## 5. Discussion

SOBHBs and IMBHs provide the opportunity to measure the effect of accretion, which might affect the GW waveform at low frequencies. Our analysis suggests that a multiband detection with LISA and a ground-based detector will be able to measure the accretion parameter of strongly accreting SOBHBs to within 50% precision for a few events. For these systems, neglecting accretion in the waveform template might lead to biases in the recovered binary parameters. These biases can be alleviated by an accurate measurement of the time of coalescence by a ground-based detector.

IMBHs in the local universe, if they exist as LISA sources, might also provide very accurate measurements of the accretion rate. Overall, for these systems, the effect of accretion should be included in the waveform to avoid bias in the intrinsic binary parameters.

Finally, accretion does not affect sky localization by LISA for SOBHBs, and it impacts that of IMBHs only mildly. In both cases, the measurement errors are typically well within Athena's and SKA fields of view. Furthermore, the X-ray flux expected from strongly accreting binaries is comparable with

Athena's sensitivity and is well above the sensitivity of future missions such as Lynx (The Lynx Team 2018). Likewise, in the case of jets, the radio signal from IMBHs could be detectable by SKA. Our analysis shows that the simultaneous operation of Athena/SKA and LISA would therefore provide the thrilling opportunity to detect the EM counterpart of highly accreting black-hole binaries.

We are indebted to Davide Gerosa, Cole Miller, Luigi Stella, and Alex Nielsen for insightful comments. A.C. acknowledges support from the “Generalitat Valencian” (Spain) through the “plan GenT” program (CIDEGENT/2018/019), as well as national grants FPA2014-57816-P, FPA2017-85985-P and the European projects H2020-MSCA-ITN-2015//674896-ELUSIVES. L.S. acknowledges that research at Perimeter Institute is supported by the Government of Canada through Industry Canada and by the Province of Ontario through the Ministry of Research and Innovation. E.B. acknowledges financial support provided under the European Union's H2020 ERC Consolidator Grant “GRavity from Astrophysical to Microscopic Scales” grant agreement No. GRAMS-815673. P.P. acknowledges financial support provided under the European Union's H2020 ERC, Starting Grant agreement No. DarkGRA-757480, and under the MIUR PRIN and FARE programmes. This project has received funding (to E.B. and A.T.) from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 690904. The authors would like to acknowledge networking support by the COST Action CA16104 and support from the Amaldi Research Center funded by the MIUR program “Dipartimento di Eccellenza” (CUP: B81118001170001).

## Appendix A: Accretion Term in the GW Waveform

An accreting binary can be described by a Hamiltonian $H(q, p)$, where the masses vary adiabatically. As shown, for instance, in Landau & Lifshitz (1960) and Sivardière (1988), the action variables $I_q = \oint pdq/(2\pi)$ are adiabatic invariants. In our case, working in polar coordinates $r, \phi$ and in the center of mass frame, we then have that $I_{\phi} = p_{\phi}$ and $I_r = \oint p_r dr/(2\pi)$ are conserved under accretion. The latter implies that circular orbits remain circular under accretion, while the former is equivalent to the conservation of the orbital angular momentum under accretion.

Then, to leading order, angular momentum is only lost through GWs (Peters 1964),

$$\dot{L}_{\text{GW}} = - \frac{32 m_1^2 m_2^2 M^{1/2}}{5 r^{7/2}} = - \frac{32}{5} \mu^2 \omega^{7/3} M^{4/3}. \quad (13)$$

Defining the reduced angular momentum $\tilde{L}_z = L_z/\mu M = \sqrt{r/M}$, the evolution of the binary can be obtained through

$$\dot{\tilde{L}}_z = \frac{\dot{L}_{\text{GW}}}{\mu M}. \quad (14)$$

Integrating Equation (14), we find the evolution of the orbital frequency,

$$\omega_{\text{GW}}(t) = \left( \omega_0^{-8/3} - \frac{256}{5} M^{2/3} \mu t \right)^{-3/8}, \quad (15)$$

<footnote n="13">¹³ The peak sensitivity is reached when the orbital velocity is equal to that of the innermost circular orbit.</footnote>
where $\omega_0 = \pi f_0$ is the initial orbital frequency. The time as a function of the orbital frequency is found inverting this expression,

$$t_{\text{GW}}(\omega) = t_c - \frac{5}{256 \mu M^{2/3} \omega^{8/3}}, \quad (16)$$

where $t_c$ is the merger time in the Newtonian approximation. In the stationary-phase approximation, the GW phase reads (Cutler & Flanagan 1994; Maggiore 2008)

$$\phi_{\text{GW}}(f) = 2\pi f t_c + \phi_c - 2 \int_{t}^{t_c} \omega_{\text{GW}}(t') dt' \\ = 2\pi f t_c + \phi_c + \frac{3}{4} (8\pi M f)^{-5/3},$$

where $\phi_c$ is the phase at merger.

We shall now compare these known results with what happens in the presence of mass accretion. We assume that the binary is surrounded by gas and that both bodies are accreting

We can now solve the total angular momentum variation equation for the orbital frequency,

$$\omega_{\text{acc}}(t) = 5^{3/8} e^{f_{\text{Edd}} \frac{(3\xi+5)}{\tau} t} \times \left( 5 \omega_0^{-8/3} - \frac{768 \mu_0 M_0^{2/3} \tau (e^{f_{\text{Edd}} \frac{(24\xi+35)}{3\tau} t} - 1)}{(24\xi + 35) f_{\text{Edd}}} \right)^{-3/8}, \quad (21)$$

where $M_0$ and $\mu_0$ are the initial values of the total and reduced mass, respectively. This expression cannot be inverted exactly to find $t = t(\omega)$. We therefore use a perturbative expansion valid when the accretion correction is small, i.e., we assume $t_{\text{acc}}(\omega) = t_{\text{GW}}(\omega) + f_{\text{Edd}} t_{\text{acc}}^{(1)}(\omega) + \mathcal{O}(f_{\text{Edd}}^2)$. We verified this to be an excellent approximation in all realistic situations, including when $f_{\text{Edd}} \sim 1-100$. This is because the dimensionless parameter always appears in the combination $f_{\text{Edd}} \cdot t/\tau$, which is always small for the evolution times scales that we consider.

In terms of the GW frequency, we find

$$t_{\text{acc}}^{(1)}(f) = - \frac{25(\pi^{16/3}f^{16/3}(24\xi + 35) - 3\pi^{16/3}f_0^{16/3}(8\xi + 15) + 10\pi^{16/3}(f_0)^{8/3})}{393216 \pi^{32/3} (f_0)^{16/3} \mu_0^2 M_0^{4/3} \tau} \quad (22)$$

mass at a same fraction of the Eddington rate,

$$m_i(t) = m_{i,0} e^{f_{\text{Edd}} t/\tau}, \quad (17)$$

where $\tau = 4.5 \times 10^7$ yr is known as the Salpeter timescale and $m_{i,0}$ is the initial mass of the $i$th body. When this time dependence is taken into account in the expression for the angular momentum, Equation (14) acquires an extra term:

$$\dot{L}_z = \frac{\dot{L}_{\text{GW}}}{\mu M} - \tilde{L}_z \frac{(\mu \dot{M})}{\mu M}. \quad (18)$$

In this equation, all masses should be considered time dependent, except the ones appearing in the angular momentum radiated by GWs. This is because accretion cannot be considered adiabatic compared to GW emission.

Accretion will in general be accompanied by a drag force $F_{\text{drag}}$ due to the fact that the accreted material carries some angular momentum. This effect can be quantified as

$$F_{\text{drag},i} = \dot{m}_i (v_{\text{gas}} - v_i). \quad (19)$$

for each mass, where $v_i$ is the velocity of the $i$th body. For simplicity, we parameterize this effect with a constant factor $\xi$, fixed by the relative velocity between the gas and the perturber (Barausse & Rezzolla 2008; Barausse et al. 2014),

$$F_{\text{drag},i} \approx -\xi \dot{m}_i v_i \quad \rightarrow \quad \dot{L}_{\text{drag}} = -\xi \dot{\mu} r^2 \omega. \quad (20)$$

Note that the parameter $\xi$ can be positive (drag) or negative (pull, see, e.g., Gruzinov et al. 2019). At leading order in $f_{\text{Edd}} \xi$, the term $\dot{L}_{\text{drag}}/\mu M$ should be added to the right-hand side of Equation (18) to take the effect of the drag into account.

Finally, we can compute the contribution of accretion to the GW phase in the stationary-phase approximation, $h \sim |h| e^{i\phi}$, at first order in perturbation theory, i.e., $\phi \simeq \phi_{\text{GW}} + \phi_{\text{acc}} = 2\pi f (t_{\text{GW}} + f_{\text{Edd}} t_{\text{acc}}^{(1)}) - \int_{0}^{f_{\text{GW}}+f_{\text{Edd}} t_{\text{acc}}^{(1)}} 2 \omega_{\text{acc}} dt$. We find, again as a function of the GW frequency,

$$\phi_{\text{acc}} = -f_{\text{Edd}} (8 \xi + 15) \frac{75 M_0}{851968 \tau} (\pi f M_0)^{-13/3} \\ + f_{\text{Edd}} \frac{25}{32768 \pi^{8/3} f_0^{8/3} M_0^{5/3} \tau} (\pi f M_0)^{-5/3} \quad (23) \\ + f_{\text{Edd}} (3 \xi + 4) \frac{25}{19968 \pi^{13/3} f_0^{13/3} M_0^{10/3} \tau} \\ - f_{\text{Edd}} (24 \xi + 35) \frac{25}{196608 \pi^{16/3} f_0^{16/3} M_0^{13/3} \tau} (\pi f M_0). \quad (24)$$

In the expression above, the terms linear in frequency and independent of frequency can be reabsorbed in the definition of the time to coalescence $t_c$ and the phase at coalescence $\phi_c$, respectively. Equation (23) tells us that the GW signal will be dominated by the effect of accretion if the frequency is sufficiently low. By comparing the size of the leading-order phase term (0PN) in the vacuum waveform and the -4PN term induced by accretion, we find that accretion is the dominant effect at frequencies below

$$f_{\text{acc}} \approx \frac{1}{\pi} \left( \frac{25}{3} \frac{45}{6656} \frac{f_{\text{Edd}}}{\tau} \right)^{3/8} M^{-5/8}. \quad (25)$$

While in Equation (23), we show all of the terms of the expansion, we have verified that the -4 PN term dominates. The inclusion of the 0PN term changes the results of the main
text by less than 1%. This is expected since most of the binary evolution in the LISA band takes place at large separation/low frequencies.

In the analysis presented in the main text, we discarded the terms proportional to the drag coefficient $\xi$, which would add an additional parameter in our waveform and require proper modeling of the distribution of the gas and its velocity around the black holes. From the functional form of Equation (23), we can see that neglecting the drag does not affect the frequency dependence of the GW phase, while it might affect the size of the effect. However, $f_{Edd}$ and $\xi$ enter the two leading terms in Equation (23) in different combinations, which would help disentangle the two effects. Indeed, we checked that for generic values of $\xi$, the time and phase shifts presented in Figure 1 do not vary dramatically.

## Appendix B
### Details on the MCMC Analysis

Using Bayes' theorem, we compute the posterior distribution for $\zeta$, the multidimensional vector parameterizing a waveform template $h$ given observed data $d$:

$$p(\zeta|d) = \frac{p(d|\zeta)p(\zeta)}{p(d)} \qquad (26)$$

For the prior $p(\zeta)$, we assume a flat distribution in $m_1$ and $m_2$ with $m_1 \ge m_2 \ge 3M_\odot$, flat in spin magnitude between -1 and 1, volume uniform for the source localization and flat in the source orientation, its polarization and its initial phase. In the LISA-only scenario, we assume a flat prior in initial frequency, and in the LISA+Earth scenario, we use instead a Gaussian prior centered around the true value of $t_c$ of width $\sigma_{t_c} = 10^{-3}$s. Assuming Gaussian noise, the likelihood is given by $p(d|\zeta) = e^{-\frac{1}{2}(d-h|d-h)}$ where parenthesis denote the inner product defined by: $(h_1|h_2) = 4 \Re e \left( \int \frac{\tilde{h}_1(f) \tilde{h}_2^*(f)}{S_n(f)} df \right)$. In the denominator, $S_n(f)$ is the detector power spectral density, indicating the level of noise at a given frequency.

To sample the posterior distribution, we use a Metropolis Hashtings Markov Chain Monte Carlo (MHMCMC; Karandikar 2006; Chib & Greenberg 1995) algorithm that we designed for this problem. More details will be given in an upcoming publication (Marsat et al. 2020; A. Toubiana et al. 2020a, in preparation). The basic idea of the algorithm is to explore the parameter space through a Markov chain generated with a symmetric proposal $\pi$, $\pi(\zeta_1, \zeta_2) = \pi(\zeta_2, \zeta_1)$. Starting from a point $\zeta_0$, we accept the proposed point $\zeta_p$ with a probability given by the ratio of the posterior distribution, $p(\zeta_p)/p(\zeta_0)$. By doing so, we accumulate samples representing the distribution. In order to increase the sampling efficiency, we parameterize the waveforms with parameters for which—based on the PN expressions (Buonanno et al. 2007, 2009)—we believe the posterior distribution is simpler. We take $\zeta = (M, \mu/M, f_0, \chi_s, \chi_a, f_{Edd}, \phi, \sin(\theta), \psi, \phi_0, \cos(\iota), \log 10 d_L)$ in the LISA-only scenario. In the LISA+Earth scenario, we use $t_c$ instead of $f_0$. Here, $\chi_s$ is the symmetric combination of spins

$$\chi_s = \frac{m_1\chi_1 + m_2\chi_2}{m_1 + m_2}, \qquad (27)$$

while $\chi_a$ is the corresponding antisymmetric combination,

$$\chi_a = \frac{m_1\chi_1 - m_2\chi_2}{m_1 + m_2}. \qquad (28)$$

For the proposal $\pi$, we use a Gaussian distribution based on Fisher matrix. To ensure we have independent samples, we downsample the chain using the autocorrelation length.

To strengthen our confidence in our MHMCMC, we cross-checked our results obtained with it to the ones obtained with Multinest, a public nested sampling algorithm (Feroz et al. 2009; Skilling 2006).

## ORCID iDs

Andrea Caputo https://orcid.org/0000-0003-1122-6606

Enrico Barausse https://orcid.org/0000-0001-6499-6263

Paolo Pani https://orcid.org/0000-0003-4443-1761

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, ApJL, 818, L22

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, ApJL, 826, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, ApJL, 882, L24

Ade, P. A. R., Agnanim, N., Arnaud, M., et al. 2016, A&A, 594, A13

Antonini, F., & Perets, H. B. 2012, ApJ, 757, 27

Atakan Gurkan, M., Freitag, M., & Rasio, F. A. 2004, ApJ, 604, 632

Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv:1702.00786

Barausse, E., Cardoso, V., & Pani, P. 2014, PhRvD, 89, 104059

Barausse, E., Cardoso, V., & Pani, P. 2015, J. Phys. Conf. Ser., 610, 012044

Barausse, E., Nicolás, Y., & Katie, C. 2016, PhRvL, 116, 241104

Barausse, E., & Rezzolla, L. 2008, PhRvD, 77, 104027

Buonanno, A., Cook, G. B., & Pretorius, F. 2007, PhRvD, 75, 124018

Buonanno, A., Iyer, B., Ochsner, E., Pan, Y., & Sathyaprakash, B. S. 2009, PhRvD, 80, 084043

Cardoso, V., & Maselli, A. 2019, arXiv:1909.05870

Carson, Z., & Yagi, K. 2020a, CQGra, 37, 02LT01

Carson, Z., & Yagi, K. 2020b, PhRvD, 101, 044047

Charles, P. A., & Coe, M. J. 2003, arXiv:astro-ph/0308020

Chen, X., & Shen, Z. 2019, MDPI Proc., 17, 4

Chib, S., & Greenberg, E. 1995, The American Statistician, 49, 327

Cutler, C. 1998, PhRvD, 57, 7089

Cutler, C., & Flanagan, É. E. 1994, PhRvD, 49, 2658

Dewdney, P., Turner, W., Braun, R., et al. 2014, SKA1 SYSTEM BASELINEV2 DESCRIPTION, SKA-TEL-SKO-0000308, https://www.skatelescope.org/wp-content/uploads/2014/03/SKA-TEL-SKO-0000308_SKA1_System_Baseline_v2-descriptionRev01-part-1-signed.pdf

Droz, S., Knapp, D. J., Poisson, E., & Owen, B. J. 1999, PhRvD, 59, 124016

Ducci, L., Sasaki, M., Haberl, F., & Pietsch, W. 2013, A&A, 553, A7

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Flanagan, E. E., & Hughes, S. A. 1998, PhRvD, 57, 4566

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics (3rd; Cambridge: Cambridge Univ. Press)

Gnocchi, G., Maselli, A., Abdelsalhin, T., Giacobbo, N., & Mapelli, M. 2019, PhRvD, 100, 064024

Gruzinov, A., Levin, Y., & Matzner, C. D. 2019, MNRAS, 492, 2755

Hartwig, T., Volonteri, M., Bromm, V., et al. 2016, MNRAS, 460, L74

Holgado, A. M., & Ricker, P. M. 2019, ApJ, 882, 39

Husa, S., Khan, S., Hannam, M., et al. 2016, PhRvD, 93, 044006

Karandikar, R. L. 2006, Sadha, 31, 81

Khan, S., Husa, S., Hannam, M., et al. 2016, PhRvD, 93, 044007

Kinugawa, T., Inayoshi, K., Hotokezaka, K., Nakauchi, D., & Nakamura, T. 2014, MNRAS, 442, 2963

Landau, L. D., & Lifshitz, E. M. 1960, Course on Theoretical Physics: Mechanics (New York: Wiley)

Latif, M. A., & Ferrara, A. 2016, PASA, 33, e051

Leigh, N. W. C., Sills, A., & Boker, T. 2013, MNRAS, 433, 1958

Lindblom, L., Owen, B. J., & Brown, D. A. 2008, PhRvD, 78, 124020

Madau, P., & Rees, M. J. 2001, ApJL, 551, L27

Maggiore, M. 2008, Gravitational Waves: Volume 1: Theory and Experiments (Oxford: Oxford Univ. Press)

Mangiagli, A., Klein, A., Sesana, A., Barausse, E., & Colpi, M. 2019, PhRvD, 99, 064056

Mapelli, M. 2016, MNRAS, 459, 3432
Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50

Marsat, S., Baker, J., & Canton, T. D. 2020, arXiv:2003.00357

Marsat, S., & Baker, J. G. 2018, arXiv:1806.10734

McGee, S., Sesana, A., & Vecchio, A. 2020, NatAs, 4, 26

McKernan, B., Saavik Ford, K. E., Bellovary, J., et al. 2018, ApJ, 866, 66

Meidinger, N. 2018, CoSka, 48, 498

Mezcua, M. 2017, IJMPD, 26, 1730021

Miller, J. M., Fabian, A. C., & Miller, M. C. 2004, ApJL, 614, L117

Miller, M. C., & Colbert, E. J. M. 2004, IJMPD, 13, 1

Miller, M. C., & Hamilton, D. P. 2002, MNRAS, 330, 232

Moesta, P., Alic, D., Rezzolla, L., Zanotti, O., & Palenzuela, C. 2012, ApJL, 749, L32

Moore, C. J., Gerosa, D., & Klein, A. 2019, MNRAS, 488, L94

Palenzuela, C., Lehner, L., & Liebling, S. L. 2010, Sci, 329, 927

Peters, P. C. 1964, PhRv, 136, B1224

Portegies, Z., Simon, F., Baumgardt, H., et al. 2004, Natur, 428, 724

Portegies, Z., Simon, F., & McMillan, S. L. W. 2002, ApJ, 576, 899

Poutanen, J., Fabrika, S., Butkevich, A. G., & Abolmasov, P. 2007, MNRAS, 377, 1187

Rodriguez, C. L., Farr, B., Raymond, V., et al. 2014, ApJ, 784, 119

Ryu, T., Tanaka, T. L., Perna, R., & Haiman, Z. 2016, MNRAS, 460, 4122

Sadowski, A. 2011, arXiv:1108.0396

Schneider, R., Ferrara, A., Natarajan, P., & Omukai, K. 2002, ApJ, 571, 30

Sesana, A. 2016, PhRvL, 116, 231102

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Sivardière, J. 1988, EJPh, 9, 150

Skilling, J. 2006, Bayesian Analysis, 1, 833

Steiner, J. F., McClintock, J. E., & Narayan, R. 2012, ApJ, 762, 104

Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946

Tamanini, N., Caprini, C., Barausse, E., et al. 2016, JCAP, 1604, 002

Tamanini, N., Klein, A., Bonvin, C., Barausse, E., & Caprini, C. 2020, PhRvD, 101, 063002

The Lynx Team 2018, arXiv:1809.09642

Vallisneri, M. 2008, PhRvD, 77, 042001

Yunes, N., & Pretorius, F. 2009, PhRvD, 80, 122003