Observing GW190521-like binary black holes and their environment with LISA

Binaries of relatively massive black holes like GW190521 have been proposed to form in dense gas environments, such as the disks of Active Galactic Nuclei (AGNs), and they might be associated with transient electromagnetic counterparts. The interactions of this putative environment with the binary could leave a significant imprint at the low gravitational wave frequencies observable with the Laser Interferometer Space Antenna (LISA). We show that LISA will be able to detect up to ten GW190521-like black hole binaries, with sky position errors $\lesssim1$ deg$^2$. Moreover, it will measure directly various effects due to the orbital motion around the supermassive black hole at the center of the AGN, especially the Doppler modulation and the Shapiro time delay. Thanks to a careful treatment of their frequency domain signal, we were able to perform the full parameter estimation of Doppler and Shapiro-modulated binaries as seen by LISA. We find that the Doppler and Shapiro effects will allow for measuring the AGN parameters (radius and inclination of the orbit around the AGN, central black hole mass) with up to percent-level precision. Properly modeling these low-frequency environmental effects is crucial to determine the binary formation history, as well as to avoid biases in the reconstruction of the source parameters and in tests of general relativity with gravitational waves.


I. INTRODUCTION
In the latest observing run, LIGO and Virgo detected a stellar-mass black hole (BH) binary, GW190521, whose features suggest the binary might have undergone stronger than expected interactions with its environment. Firstly, GW190521 had component masses 1 [3,4], with the larger lying squarely within the pairinstability gap ∼ [50,130]M [5][6][7]. This has prompted a lively debate in the literature, with suggestions that the progenitor BHs may have formed not via standard stellar evolution in the field, but rather by repeated coalescences in dense environments. The latter include globular or nuclear stellar clusters [8-10] or active galactic nuclei (AGNs) [11][12][13], where BHs may migrate to and accumulate in the nuclear region at faster pace than stars due to their larger mass. This in turn enhances their growth by mergers and accretion. These "dynamical" formation channels would also explain the apparently large and misaligned spins of the component BHs of GW190521 [3,4].
Secondly, the Zwicky Transient Facility (ZTF) detected an optical flare (dubbed ZTF19abanrhr) about 34 days af-ter GW190521, in AGN J124942.3+344929 at redshiftz s = 0.438. The position and distance of this system are compatible with the inferred position and distance of GW190521, especially under the assumption of a uniform prior in massratio [2]. Reference [14] then interpreted ZTF19abanrhr as due to the BH remnant from GW190521 moving in the AGN disk (as a result of the recoil produced by the anisotropic gravitational wave emission during the merger). In this picture, the 34-day delay between GW190521 and ZTF19abanrhr would be due to the time required by the radiation produced at the shock front induced by the recoiling BH remnant to emerge from the surface of the disk. The recoil has been estimated at a velocity ∼ 200 km/s at ∼ 60 deg from the midplane of the disk, whose aspect ratio -height H to galactocentric radius a • -was inferred to be H/a • ∼ 0.01. Reference [14] also argued that the GW190521 binary was likely located in a disk migration trap (where gas torques cancel out and binaries accumulate in their inward migration [15]), i.e., its distance from the nucleus should be a • ∼ 700GM • /c 2 , with M • ∼ 10 8 − 10 9 M the mass of the BH at the center of the AGN.
If indeed GW190521 lived in an AGN disk, it could be the harbinger of a significant population of binary BHs dwelling in gas-rich environments that will be detected in the coming years by ground and space detectors. The large gas densities and the presence of supermassive BHs (SMBHs) at the center of AGN disks could produce observable effects in gravitational wave (GW) signals. For instance, while the gravi-tational pull of an AGN disk is negligible [16], the gravitational interaction of a BH with the sound waves (wakes) that it excites in the surrounding gas leads to effects (dynamical friction and planetary-like migration) that can affect the longterm evolution of a BH binary system [17][18][19][20]. The surrounding gas can also accrete onto the component BHs, transferring linear momentum to them [17,19,20]. Moreover, the presence of the central SMBH can significantly affect the waveform, as a result of the accelerated motion of the stellar-origin BH binary around it [21][22][23], or because of lensing [24,25] and Shapiro time delay [26].
The most significant among these environmental effectsaccretion, dynamical friction, acceleration -are believed to be especially prominent at low frequencies, as they appear at negative post-Newtonian (PN) orders in the GW phase relative to the vacuum quadrupole emission [19,27]. While current ground based detectors experience a degradation of their sensitivity at low frequencies (due to seismic noise), environmental effects could be in principle observable with the Laser Interferometer Space Antenna (LISA), by targeting the early inspiral of stellar-origin BH binaries months/years before they merge in the band of ground detectors [20,22,23,[27][28][29][30][31][32]. This opens the possibility to use the detector to shed light on the physics of AGN disks and on the properties of SMBHs at their centers [31,32], as well as to gather important clues on the formation mechanisms of stellar-mass BH binaries [33].
On the downside, matter effects could either be degenerate with tests of general relativity (GR) with low-frequency GW signals, or bias the results of such tests if not properly taken into account. Indeed, some modifications of general relativity also predict low-frequency phase contributions, some already partially constrained by LIGO/Virgo observations [34,35]. LISA has the potential to improve low-frequency constraints by about 8 orders of magnitude (in the flux) [36,37], provided that matter effects are absent or accurately modelled.
In this work, we extend our previous analysis of GW190521-like binaries [27], focusing on the modifications to the vacuum waveform induced by the third body (the SMBH) and assessing the detectability of its effect with LISA. While motivated by GW190521 and the AGN formation scenario, our results apply more broadly to stellar-mass binaries orbiting a third body with a period comparable to the GW observation time, in the absence of gas.
Three body systems have a rich and varied phenomenology, depending on the hierarchy of masses and distances of the three bodies. Here, we focus on hierarchical triples in which a compact, circular stellar-mass binary is on a circular orbit around a SMBH with period comparable to LISA's observation time. For these systems, Refs. [24,28,31,38] found that the Doppler modulation could be measured and used to infer the SMBH mass. If the outer orbit is sufficiently aligned with the line of sight, Refs. [24,39] showed that lensing could also be detected by LISA. If the binary is even closer to the SMBH, its spin can also be measured thanks to the Lens-Thirring effect [31], and retro-lensing [39] can be significant. Finally, Refs. [29,30,40] explored other, higher order effects (such as the Kozai-Lidov eccentricity oscillations) and found that these could also be detected by LISA.
These studies made a number of simplifying assumptions, in modelling the Doppler-modulated waveform in the frequency domain with the stationary phase approximation (SPA) and/or by relying on the (sometimes approximate) Fisher matrix to estimate LISA's constraining power, and neglecting the Shapiro time delay. In this work, we address these limitations, and perform realistic parameter estimation for the first time for a binary orbiting a SMBH. We provide a recipe to obtain the frequency domain waveform for signals strongly modulated by the peculiar motion and the Shapiro delay, improving on the global SPA commonly used in the literature. With this improved and fast model, we perform a full Bayesian analysis and show that LISA will be able to measure the SMBH mass with percent precision in GW190521like systems.
This paper is organized as follows: in Section II we briefly summarize how many (vacuum) GW190521-like binaries LISA could detect, and how well it could measure their vacuum parameters. In Section III we describe how the waveform is modified by the two dominant third-body phenomena: the Doppler and Shapiro effects. We also review other thirdbody effects, finding them negligible for the system configuration suggested by the ZTF counterpart. Finally, in Section IV we model the Doppler and Shapiro-modulated signal in frequency domain and present the results of parameter estimation. We also include some results on the detectability of these signals with vacuum templates. We conclude with some prospects in Section V.
Throughout our work, we use the following notation: m 1 , m 2 are the masses of the inner (GW190521-like) binary,z s the cosmological redshift of the source, M • the mass of the SMBH at the centre of the AGN, ι • the inclination of the outer orbit with respect to the line of sight, a • the radius of the outer orbit, f the frequency of the GW. For the outer orbit parameters, see also Fig. 2. We work in units in which G = c = 1.

II. OBSERVING (VACUUM) GW190521-LIKE BINARIES WITH LISA
In this Section, we briefly summarize the prospects of detecting and inferring the properties of GW190521-like binaries with LISA, assuming their vacuum GW signal.

A. Rates and detectability in isolation
We compute the average number of expected events according tō (2.1) where θ stands for the angles (inclination, position in the sky and polarization), t c is the observed coalescence time, z s is the source redshift, m 1 > m 2 the source-frame component masses, V c the comoving volume, R the intrinsic rate (assumed constant in redshift, details below) and p(m 1 , m 2 ) noise, T obs , DC GW190521-like GWTC-3 GWTC-3 massive  the mass probability density. We take angles, merger time and comoving volume to be uniformly distributed and assume the Planck 2015 ΛCDM cosmology [41]. We explore different detection thresholds for the signal-to-noise ratio SNR thr , mission duration times, instrument duty cycles (SNR DC = SNR DC=100% × √ DC) and instrument noise curves: "SciRD", with OMS (optical metrology system) noise of 15 pm [42] and "MRD", with OMS noise of 10 pm [43] 2 . In computing the SNR, we neglect the motion of the instrument. We focus on events betweenz min = 0 andz max = 0.512 (corresponding to a comoving distance of 2 Gpc in our cosmology) and merging within 100 years from the beginning of the LISA mission, as environmental effects can be measured better in events with shorter merger times [20].
We explore two different ways of predicting the intrinsic rate. In the first case ("GW190521-like" in Table I) we draw realizations of masses from the posterior distribution obtained by the LIGO/Virgo collaboration (LVC) when analyzing GW190521 with the NRSur waveform model [4]. We draw the rate from a Gamma distribution whose median and 90% CI approximately match the values 0.08 +0. 19 −0.07 Gpc −3 yr −1 , the latest LVC estimate for GW190521-like mergers assuming no redshift evolution [44]. For each rate R, we produce samples of the number of events from a Poisson distribution with meanN computed according to Eq. (2.1). Finally, we compute the median and symmetric 90% confidence interval after sampling over both the rate and the universe realization. The results are shown in Table I, and are in agreement with Ref. [45]. We conclude that LISA could detect a few events like GW190521, with a large uncertainty on the actual number depending on the poorly constrained rates and on the actual LISA mission configuration.
In the second case ("GWTC-3", in Table I), we use the rate estimate (with median 28.1 Gpc −3 yr −1 ) obtained by the LIGO/Virgo/KAGRA collaboration by analyzing the O1, O2, O3a-b events, assuming a uniform merger rate in comoving volume-time and a "power law+peak" population model [46]. We sample from the posterior distribution of the power law+peak model parameters provided in Ref. [46]. For each population model realization we estimate the mean number of events, and produce realizations from a Poisson distribution with meanN. Finally, we compute the median and symmetric 90% confidence interval and report them in Table I.
The formation mechanism of the GWTC-3 events is still being debated, with one Bayesian analysis finding for about 25% of them a preference for an AGN origin (in a specific model) [47]. To mimic the AGN sub-population, we consider the sub-set of events with primary masses above 45 M (the column named "GWTC-3 massive" in Table I).
Our predictions for the overall rate of stellar-mass binaries detectable by LISA lie between a few and O(100), and are broadly consistent with other estimates in the literature once we account for differences in the rate, detection threshold, etc. (see, e.g., [48,49]). 3 Furthermore, based on the latest catalog, we predict there will be O(10) events detected by the LISA mission likely to occur in AGNs (still with a considerable uncertainty). Note that, within the detectable events in any of the configurations in Table I, only a fraction will be realistic (merging within 10 years) multiband events.
In Fig. 1 we also show the distribution of the 10-year SNR of GW190521-like binaries in LISA, averaged over the sky position and polarization angle, and with merger time of 10 yrs. Figure 1 also reports the number of GW190521-like events detectable by LISA as a function of the SNR threshold, for a 75% duty-cycle and a fiducial (SciRD, 6 yrs) and optimistic (MRD, 10 yrs) mission configuration.
Overall, the SNR distribution in Fig. 1, as well as the rates in Table I, point to the fact that GW190521-like binaries will be hard to detect with LISA, unless found at a closer distance. Detection prospects would improve with strategies (such as multiband detections [48]) requiring a lower SNR threshold, as seen in Fig. 1, right panel.

B. Parameter estimation in isolation
Despite the large masses of GW190521-like systems, (as defined by the samples released by the LVC), the morphology of these GW signals is very similar to the one of stellarmass BH binaries, i.e., we observe the system during its early inspiral and it then leaves the LISA band as it starts chirping. We therefore refer to Ref. [32] for a complete description of the correlations between parameters and the accuracy in their measurement. Among intrinsic parameters only the chirp mass can be precisely measured, with relative error 10 −4 . Thanks to its long duration in the LISA band, the source can be very well localised, within 1 deg 2 , see also Ref. [45]. Finally, the distance to the source is measured with a 60% error. As we will see in Section IV, the measurement of the binary parameters barely changes when it is near a SMBH.
FIG. 2. Snapshot of the system when the inner binary is behind (at t = t 1 ) or in front (at t = t 2 ) of the SMBH, i.e., when the source, the central black hole and the observer all lie in the page plane. The outer orbit is perpendicular to this plane.

III. GW190521 AROUND AN AGN
In Ref. [27] some of the current authors assessed the detectability of deviations from purely vacuum waveforms arising from AGN matter effects, such as gas accretion and dynamical friction. If GW190521-like events occur in dense gaseous environments, with sizable values for the gas density, ρ gas 10 −10 g/cm 3 , and accretion rates, f Edd O(1), then LISA will be able to detect these effects in the GW signal. In the following, we will focus our attention on the AGN environmental effects related to the presence of a third body, which have been not analyzed in depth in previous works.
We will focus on binaries at 1) intermediate distances from the SMBH, with a • = O(100)M • , and 2) with generic inclinations between the outer orbit and the line of sight. Because of the putative identification with an electromagnetic counterpart [14], we also refer to these binaries as GW190521-like.
The two dominant effects in these conditions are the Doppler and Shapiro effects, which we describe in detail in this Section. We will only briefly discuss effects that might be relevant at shorter separations (e.g. Lense-Thirring precession) or for highly aligned systems (lensing).
Our choice for the distance from the SMBH is motivated by the fact that we expect AGN binaries to be preferentially located in disk traps [51,52], where inward and outward disk torques balance and black holes can accumulate and merge hierarchically. These traps are typically found between O(10)M • and O(10 3 )M • from a central SMBH of mass M • [53] (although, see also [54]), and could even occur close to the innermost stable circular orbit [55].
We also assume that the outer orbit is circular, as migrating bodies in AGN disks, similarly to planets in protoplanetary disks, are expected to circularize [56,57]. Finally, at distances O(100)M • from the central object, we can assume that the SMBH is nonrotating (spin-related effects are negligible, see Section III D 2).
For the inner binary, we restrict our analysis to a spinaligned and circularized system, leaving the modelling of spin-induced precession and eccentricity in binaries around a SMBH (see, e.g., [58,59]) for future studies.

A. The waveform in the observer frame
We now derive the waveform in the observer frame ("o") and relate it to the usual expression in the source frame ("s"). The relation between the observer-frame and source-frame time, as well as the observed redshift, can be derived as follows. We use the McVittie metric [60,61] with µ(t, r) = m/(2ra(t)) to describe an observer in an expanding (flat) Universe in the presence of a static object of mass m (corresponding to a Schwarzschild BH when a(t) =1), from which the coordinate r originates. This is a good description for sources at cosmological distances ( Gpc) from us -for those at closer distances, we could simply use the Schwarzschild metric. For m/r 1, Eq. (3.1) can be put in the standard form of a perturbed expanding universe in longitudinal gauge (see, e.g., [62]), and we neglect the motion and the potential at the observer, which can be taken into account separately, in the GW detector response. At the source, the wave-vector of the gravitational wave, perturbed by the presence of the point mass (i.e. the central BH), is k µ s = a −2 (1 + δk 0 s ,n + δ k s ) withn pointing from the source to the observer, and where we again neglect perturbations at the observer. Because we assumed a circular outer orbit, the magnitude of the source velocity is constant in time, γ s = const. The redshift between the observer and the source is then given by [62,63] ( where we neglected higher order terms mixing the potential and the peculiar motion, as well as the time derivatives of the potentials. We also absorbed the constant γ s factor in the cosmological redshift γ s a o /a s → (1 +z s ). From Eq. (3.3), we see that the observed redshift receives contributions from the cosmological redshiftz s , the peculiar motion (Doppler effect, z D ) and the gravitational influence of the third body (Shapiro effect, z S ). We can integrate the redshift to obtain the time measured by the observer. We neglect the change in the expansion of the universe during the time the signal is propagating from the source to the observer, so that 1 +z s = const (for works taking into account the variation of the background expansion, see [21,23,64,65]). Integrating from a reference time that we set to zero for simplicity, we find where we define dt s,z ≡ (1+z s ) dt s . In the last line of Eq. (3.4), we expressed the delays in terms of t s,z , which takes into account the background cosmology. This is a convenient step for computing the waveform in the observer frame, as we will use waveform models that already take into account the standard propagation in a cosmological background. The second term in Eq. (3.4) represents the delay due to the peculiar motion (Doppler delay, in short, Section III B), while the third term is the gravitational delay due to the central mass (Shapiro delay, Section III C). Finally, the GW signal emitted by the binary in orbit around the SMBH and seen by our detector can be written as . Neglecting amplitude corrections, the observed signal is related to the signal in the source frame by a transformation of the argument, Here A s and ϕ s are the amplitude and phase of the GW signal taking into account the propagation from the source to the observer in the cosmological background (i.e., they are expressed in terms of the luminosity distance and redshifted masses). The observed time t o is given by Eq. (3.4).
In this work, we neglect the amplitude corrections the GW receives in the observer frame from the peculiar motion and the potential at the source, which are of the order (1 + z D,S (t s )) [24,66]. In other words, we only account for the transformation in the argument of the waveform (i.e., the time shift). We justify neglecting amplitude corrections in the sections to follow. We also focus on the l = m = 2 mode, and neglect the mode mixing generated by peculiar motion [66][67][68].
In the remainder, we describe in more detail the two effects contributing to the time shift (3.4). For simplicity, we drop the subscript in the source-frame time, setting t s,z = t.

B. The Doppler effect
The leading order effect of the outer orbit is a variation of the distance between the detector and the source, causing a varying Doppler delay d D (t).
The signal in the observer frame is given by Eqs. (3.5), (3.6), with t o = t + d D (t). The delay d D (t) is the time dependent change in the arrival time of the signal compared to a source with negligible proper motion: from the second term in Eq. (3.4), one gets (recalling that t stands now for t s,z ) where the distance r(t) is the projection onto the line of sight (pointing now from the observer to the source) of the outer orbit, which can be written as   .7), we see that the Doppler delay has a mild dependence on the inclination of the outer orbit, ι • . We also see that the distance from the SMBH enters the waveform in combination with the orbital inclination, in the factor a • cos ι • . The SMBH mass also only enters in combination with the orbital separation, in the orbital frequency Ω. When the Doppler effect is the only orbital effect measured in the GW signal, these two combinations of the orbital parameters are the only observable quantities.
At a typical disk trap location, the observed outer orbital period is which is comparable with the LISA mission duration, and we expect the delay to induce appreciable variations in the GW phase. An estimate of the amplitude of the phase modulation induced by the Doppler effect is given by This phase modulation is, at realistic orbital separations, the largest non-vacuum effect in AGN binaries, much larger than the phase modulation of the LISA spacecraft motion itself (approximately 30 rad). The large dephasing induced by the peculiar motion in time domain can be seen in Fig. 3. 4 The detectability of the peculiar motion phase shift was explored, for more moderate accelerations, in several studies [21,23,28,38].
As anticipated in Ref. [27], for orbital periods comparable or smaller than the LISA observational time, the Doppler effect modulates the signal to the point that the GW frequency evolution is not monotonic during observations, see Figs. 3 and 4. Because the time-to-frequency map is multivalued, a global use of the SPA (appearing e.g. in [31,39]) would not be adequate to describe the signal in the frequency domain. In this work, we will focus on this regime and explain how to obtain a more accurate frequency-domain waveform (Section IV).
The Doppler effect will also modulate the GW amplitude through a time dependent redshift factor. From Eq. (3.3), and keeping track of the redshift factors in the transformation from the source-frame time, This will induce a modulation of the GW amplitude similar to the modulation of the GW frequency, with amplitude δz D = z D,max − z D,min of order (3.12) Gravitational wave detectors are more sensitive to phase, rather than amplitude, modulations. For reference, in the case of a GW190521-like event at 1.5 Gpc distance, LISA will measure the luminosity distance with precision of order 50%, see Section IV, Fig. 6. A Doppler-induced amplitude modulation of percent level will thus be negligible and in any case subdominant compared to the phase modulation induced by the same effect, as we confirm numerically in Section IV. We therefore neglect the amplitude modulation when performing parameter estimation. For sufficiently large outer orbital separations, i.e., for orbital periods much longer than the LISA observational time, the peculiar motion reduces to a constant peculiar velocity and a constant centripetal acceleration, see Appendix A. In this case, the constant peculiar velocity can be reabsorbed in the constant redshift (and the latter in the redshifted chirp mass) and the centripetal acceleration produces a single -4PN term in the GW phase. Parameter estimation in this limit was discussed in Ref. [27] (see also [23]).

C. The Shapiro effect
The second most relevant effect in our setting is the Shapiro effect, due to the non-vanishing gravitational potential of the third body.
The Shapiro time delay is the delay in the time of arrival of a signal (with respect to the signal traveling in flat spacetime) due to the non vanishing gravitational potential along the line of sight. Assuming the central object is described by the Schwarzschild metric 5 , this is given by [69,70]  where κ = ±1 depending on whether the source is behind or in front of the SMBH, respectively. Here r s = 2M • is the Schwarzschild radius of the SMBH, R is the distance between the observer and the SMBH and is the impact parameter of the unperturbed trajectory from the source to the observer. See again Fig. 2 to visualize the geometry of the system. In the second line of Eq. (3.13), we used the fact that R b and that its constant contribution can be re-absorbed in the definition of the time of coalescence.
Equation (3.13) shows that the Shapiro delay introduces a direct dependence on the SMBH mass M •,z s , through r s . This breaks the degeneracy introduced with the Doppler delay (3.7) between the outer orbit inclination and the central mass. A detectable Shapiro delay will then allow LISA to measure all three environmental parameters (the SMBH mass M •,z s , inclination ι • and distance from the SMBH a •,z s ), as we show in the parameter estimation example in Sec. IV. The Shapiro delay in our system happens to be numerically similar to the light crossing time of the Earth orbit, so its modulation of the GW phase is comparable to the Doppler modulation induced by the LISA spacecraft motion, The Shapiro effect is therefore subdominant compared to the Doppler effect in GW190521-like binaries in LISA, but still detectable, see also Fig. 3. As we will see in Section IV, the Shapiro effect is strong enough to break the degenaracies between the outer orbital parameters and improve our ability to constrain the SMBH mass. From Eq. (3.15), we also see that the Shapiro phase modulation has a stronger dependence on the inclination of the outer orbit ι • compared to the Doppler effect, with the effect peaking when the orbit and the line of sight are aligned [71]. This can be seen comparing the two panels in Fig. 4. For high alignments, however, strong lensing is not negligible, as we argue below.
For completeness, we also write the Shapiro redshift by differentiating Eq. (3.13) and keeping track of all redshift factors, .

D. Other effects
In this work, we focus on systems dominated by the Doppler and Shapiro effects: GW190521-like binaries on a circular orbit of hundreds of Schwarzschild radii around a 10 8 -10 9 M SMBH. In general, hierarchical triple systems can be affected by a number of other dynamical effects, which might leave a detectable imprint on the GW signal for more extreme system parameters [24,28,31,40,67]. Here we briefly summarize these effects.

Lensing
We treat the SMBH as a point-like lens [24], with the binary moving on a sphere of radius a • centered on the lens, see Fig. 2. Projected along the line of sight, the angular-diameter distance between the observer and the SMBH and the observer and the source are denoted by D L and D S , respectively. The The second derivative of the GW phaseω ≡ d dt ω = d 2 dt 2 ϕ(t +d effect (t)) in isolation and in the presence of the Doppler and Shapiro effects. We compare a system with large (ι • = π/6 rad = 30 deg, left) and small (ι • = −0.05 rad −3 deg, right) outer orbit inclination. The remaining system parameters and plot features are the same as in Fig. 3. The figure on the right includes a zoom-in of the region marked by vertical lines, corresponding to a month around the moment the source is aligned behind the SMBH. At each passage, the Shapiro effect is strongly enhanced in the small inclination case. Note that at small inclinations the effect of lensing, not included here, would also be significant. angular-diameter distance between the lens and the source, D LS , can take either sign during the evolution of the binary on the outer orbit, and coincides with r(t) defined in the previous Section. In this lensing system, the Einstein radius is given by (assuming D S ∼ D L as a • D L , and D LS > 0), and its maximum is at Ωt One can identify three lensing regimes based upon the motion of the source with respect to the lens over its observable time [24]. The three relevant scales are the orbital period of the outer binary T , the time in band T obs and the time for the binary source to cross the Einstein radius of the lens, (3.20) The resulting three lensing regimes are 1) the repeating lens regime, for T < T obs ; 2) the slowly-moving lens regime, T ≥ T obs ≥ T lens ; 3) the stationary lens regime, T ≥ T obs ≤ T lens . GW190521-like binaries with outer orbital period given by Eq. (3.9) fall in the repeating lens or in the slowly-moving lens regime, depending on T obs and M •z s . A significant lensing event occurs when the source passes within one Einstein radius of the lens. Therefore, the probability of observing a lensing event after observing the binary for an entire orbit is the ratio of the binary inclination angle for which the source falls within one Einstein radius of the lens to the total possible range of inclinations. This probability depends on three parameters: T obs , a • and M • . In Ref. [27], we found that the probability of lensing occurring during a LISA observation was above 10% in a good portion of the lens parameter space, e.g., for a • 200M • for M • = 10 8 . However, in the remainder of this work, we decided to focus for simplicity on the most likely scenario, in which lensing (including "retro-lensing" [39]) is suppressed. This occurs for most outer orbital inclinations ι • , for which the source is always outside the Einstein radius of the lens:

Relativistic and tidal effects
Other third-body effects include: the outer orbit eccentricity, the geodetic or de Sitter precession of the inner orbit around the orbital angular momentum of the outer binary; Kozai-Lidov and other resonant oscillations of the inner binary inclination and eccentricity; the Lense-Thirring precession of the inner orbit around the SMBH spin, and aberration. For a summary of potentially relevant time delays, see e.g. [72].
To estimate the significance of these effects, one can compare their period to the observation time. While the Doppler and lensing/Shapiro effects oscillate with the orbital period, Eq. (3.9), the de Sitter precession period is [31] T dS = 8 ×   where we introduced the outer orbit eccentricity e and the SMBH spin χ • . The Kozai-Lidov oscillations depend also on the GW frequency [31], with a typical timescale given by (3.24) The estimates above show that higher order effects can be safely neglected in our systems, as they occur on a much longer timescale -unless the eccentricity of the outer orbit is high, which is not expected for binaries migrating in AGN disks. Aberration [67,73] is also suppressed by the ratio between the outer orbital period and the inner one (for the phase shift) and by the ratio between the orbital velocity and the speed of light (for the amplitude correction). Apsidal precession resonance [58], which can drive eccentricity oscillations in the inner binary, are also negligible when the outer orbit is circular (as assumed in this work in view of the circularizing effect of the AGN environment).

IV. PARAMETER ESTIMATION IN AGN
We now investigate the impact of a central SMBH on the GW emission from GW190521-like binaries detected by LISA, using the waveform model described in the previous Section. We consider a GW190521-like binary located close to the SMBH (as suggested by Ref. [14]), namely at a distance a • = 700M • for an observed SMBH mass M • = 10 8 M (see also Table II). In this Section, we drop the redshift subscript from all quantities, for simplicity. We set the inital phase of the outer orbit to φ • = −3π/4, so as not to produce a fine-tuned enhancement or suppression of the environmental effects in the limit of large orbital separation. In this exercise, we neglect the effect of the AGN disk gas and concentrate only on the Doppler modulation of the GW signal due to the orbital motion of the binary around the SMBH, as well as the timedependent Shapiro delay. We choose a generic inclination of the orbit with respect to the line of sight, so that lensing is negligible. We also neglect other possible interactions of the binary with the SMBH that are suppressed at intermediate separations (Kozai-Lidov resonances, orbital precession, etc.).

A. Waveform model and data analysis methods
As the source-frame GW model, we take PhenomD [74,75], which describes the inspiral-merger-ringdown signal for 104.28 1.68 -0.23 -0.23 1. 4   TABLE III. The (cosmologically redshifted) parameters of the GW190521-like binary selected for parameter estimation, consistent with the posterior distribution of GW190521 as analyzed by the LVC [4] and with SNR= 9.5 in LISA. The other signal parameters are ι = 0.85 rad, λ = 5.668 rad, β = −0.15 rad, ψ = 1.2 rad, φ obs . The initial GW frequency is f 0 = 0.0061524060, corresponding to t c = 7 yr in the vacuum injection. The corresponding cosmological redshift is z s 0.27 [41].
spin-aligned circularized binaries using only the dominant GW mode. For LISA observations, only the inspiral part of the signal is relevant, and precession and sub-dominant modes are severely suppressed by the post-Newtonian factor. The waveform is generated directly in the frequency domain.
As described in details in Section III, the phase of the GW will be strongly modulated by the Doppler and Shapiro effects, so that the signal will appear as chirping (when the binary is behind the SMBH) and anti-chirping (when the binary is in front of the SMBH). The orbital period of our outer binary is about 1.8 years, and assuming 6 years of LISA observations and that the binary is observed 7 years before the merger, we expect to see approximately 3 full orbital cycles. One can clearly see the three orbital cycles and two (anti)chirping phases in the signal in Fig. 4.
We transform the source-frame signal into the time domain, then model the phase modulation as a time-dependent delay, given by Eqs. (3.4), (3.7) and (3.13). We then split the signal into chirping and anti-chirping parts and transform the signal back into frequency domain using the SPA. Our piecewise SPA is an improvement over the global approximation adopted, e.g., in Refs. [31,39]. Our approximation only breaks at the turning points, which would require a higher order SPA [76]. We decided to neglect these parts of the signal as they contribute only a very small fraction of the total SNR.
In Section III, we also argued that the amplitude of the GW will receive a correction and oscillate with time as a result of the peculiar motion. This amplitude correction has a more modest effect on the waveform compared to the phase modulation: the overlap between a signal with and without amplitude correction is 0.9998 for our system parameters (and a signal with SNR= 9.3). Neglecting the amplitude modulation should introduce a bias of approximately 2% in the luminosity distance, much smaller than the measurement error for this parameter. For this reason, we neglect the amplitude correction and only model the time delay.
Before we proceed to parameter estimation, we introduce the signal parametrization. For the vacuum waveform, we use the chirp mass M c , the mass ratio q = m 1 /m 2 ≥ 1, and the following combination of physicals spins χ PN = η (113q + 75) χ 1 + (113/q + 75) χ 2 /113, χ − = q/(1 + q)χ 1 − 1/(1 + q)χ 2 , where η = q/(1 + q) 2 is the symmetric mass ratio and χ 1,2 are spins of two BHs. The other vacuum parameters are standard: the luminosity distance d L , the sky position as ecliptic longitude and latitude (λ, β), inclination ι, the polarization ψ and the azimuthal position of the observer φ obs . An important characteristic of the system is the time to merger, which we parametrize through the initial GW frequency f 0 , which is the observed GW frequency (uniquely defined in the absence of higher order modes and precession) at the start of LISA observations, and kept the same when we place the binary around the SMBH or in vacuum. The AGN-related parameters we use to parametrize the signal are the orbital (Keplerian) angular velocity Ω, the projections of the orbital radius a • cos ι • , a • sin ι • and initial phase φ • . For the GW190521-like binary that we use as a representative of the AGN binary population, we choose the parameters given in Table III. This particular choice of the parameters, over all the possible values given by the LVC posterior samples for GW190521 [4], was made to maximize the vacuum SNR. For the same reason, we also set t c = 7 years.
We perform parameter estimation using parallel tempering Markov-chain Monte-Carlo within a Bayesian framework, see [32,77] for details. The signal we have chosen produces SNR= 9.5 in vacuum, (SNR= 9.3 in the presence of the SMBH) in the LISA band using the noise budget outlined in the LISA science requirement document, SciRD [42].

B. Measurement of the central black hole orbit and properties
The complete results of the parameter estimation are given in Appendix B. Here we focus on the most interesting system parameters, marginalising over the others.
The combination of a very strong Doppler modulation and a time-varying Shapiro delay allows us to constrain the orbital parameters of the outer binary and the mass of the AGN BH, as shown in Fig. 5. For this particular system, we can determine the mass of SMBH to about 8% level and the orbital inclination within a few degrees, together with the separation of the outer orbit within 3%. Comparably precise SMBH mass measurements are currently possible only for SgrA * [78], M87 [79], and an handful of galaxies with detected nuclear megamaser emission [80]. GW190521-like binaries detected by LISA would therefore provide a competitive, complementary opportunity to measure SMBH masses through GW observations. Electromagnetic-and GW-based measurements would be affected by different systematic uncertainties, which could help strengthen the measurement accuracy.
The Doppler modulation of the phase of the GW signal is a strong effect, but it is mostly orthogonal to the vacuum phase evolution; as a result, the parameters of the AGN orbit are uncorrelated to the parameters of the emitting binary. To show this explicitly we overplot in Fig. 6, for the vacuum parameters, the posteriors obtained in the presence and in the absence of the SMBH. The only parameter for which the posterior distribution is affected by the presence of the SMBH is the initial GW frequency: the difference between vacuum (green contours) and the binary orbiting the SMBH (grey contours) is clearly seen in the figure. Note that the initial frequency is the observed one and, therefore, the initial orbital frequencies of the inner binary in the two cases are quite different, due to the Doppler shift caused by the orbital motion around the SMBH. The two binaries therefore merge at different times, although the shape of the posterior for the time of coalescence is not strongly modified by the presence of the SMBH.

C. Wider orbits and detectability
As we move the binary further away from the SMBH by increasing the size of the orbit, we approach the regime where the acceleration projected on the line of sight is constant, presented in our previous publication [27]. In this limit, the acceleration enters the GW phase at -4PN, a term otherwise absent in vacuum templates. For even broader orbits, we reach the regime where the acceleration is negligible and velocity is (almost) constant. The latter is (almost) degenerate with a small modification in the redshift of the source, causing a slight shift in the observed (redshifted) masses of the binary. Therefore, for sufficiently wide orbits, we expect a vacuum template to be a good approximation of the signal, and a good template to detect the source with matched-filtering.
To investigate the detectability of such a system using vacuum templates, we compute the fitting factor (i.e., the overlap maximized over all parameters). The fitting factor (FF), or rather (1 − FF), gives the fractional loss in SNR due to mismatch between the best-fitting model and the signal [81]. The fitting factor takes into account that we still might match the signal at the expense of a bias in the parameter space, which  can (at least partially) compensate for the mismatch in the model. We start with the binary orbiting the SMBH on a very broad orbit (a • 1 pc), where the orbital speed is almost constant and can be re-absorbed into the redshift, with FF close to 100%. In Fig. 7, the top left panel shows the evolution of the FF as we bring the binary closer to the SMBH. The top right plot gives the value of the chirp mass which maximizes the overlap between the signal and the model. In particular, we show the difference between the best-model chirp mass and either 1) the true (cosmologically redshifted) chirp mass or 2) the true chirp mass corrected by the Doppler redshift computed at the start time of observations. The lower row of plots shows the best-model mass ratio (bottom left) and the difference between the best-model time of coalescence t c and true value in the observer frame (bottom right).
The bias in the parameters recovered by our search has a clear increasing trend as we decrease the separation from the SMBH. At the closest distance we explored, a • 0.35 pc, the recovered parameters (most noticeably, the mass ratio) display a jump. This is most likely due to non-optimal recovery of the FF maximum. Below this distance, the search for the FF max- FIG. 7. Analysis of a GW190521-like binary at large distance from the SMBH with a vacuum template. We show the fitting factor (upper left), the difference between best-fit and true (cosmologically redshifted) chirp mass (upper right), best-fit mass ratio (lower left) and difference between the best-fit and true (observed) time to coalescence (lower right), as a function of the distance of the GW190521-like binary from the SMBH. For the chirp mass, we also show the difference between the best-fit value and the true value redshifted by the Doppler effect, as measured at the start of observations. As the distance from the SMBH decreases, the vacuum template causes a significant loss in SNR and strong biases in the inner binary parameters. The true value for the chirp mass and mass ratio is 104.28 M and 1.68, respectively, the measurement precision (half of the 90% confidence interval in Fig. 8) is approximately 5 × 10 −3 M for the chirp mass and 20 minutes for the time to coalescence. imum (in an 11-dimensional parameter space, although not all parameters are equally important) becomes challenging.
We find that the bias in the chirp mass, although small in an absolute sense, is quite large compared to the typical measurement error for this parameter (10 −3 M , see Fig. 6). At these large separations, however, most of the bias can be explained in terms of a constant Doppler redshift, as can be seen comparing the two curves in the top right panel of Fig. 7.
The bias in the coalescence time is also large compared to the precision with which this parameter is measured -O(20) minutes, worse than for typical stellar-origin BH binaries in LISA [32]. If such a binary was detected and analysed in LISA with vacuum templates, the bias in the coalescence time could prevent the association with the merger part of the signal detectable by ground-based detectors (unless other source parameters, such as chirp mass and sky position, already establish an association). A bias in the merger time would also affect tests of GR, as emission into non-GR polarizations can also contribute to a shift in t c (see Ref. [36]) and thus be degenerate with the environmental effect studied here.
We should stress that large biases might affect only a small fraction of AGN binaries, whose separation from the SMBH is sufficiently large for vacuum-template detection, but small enough to induce a sizeable bias. The severity of this problem for LISA will depend on how binaries are distributed within AGNs, and on the details of multiband strategies implemented by both space-borne and ground-based detectors in the future.
We also expect the luminosity distance of the source to be significantly biased [23] away from its background, cosmological value, with part of the bias explained by a constant Doppler redshift. However, since we neglected subleading amplitude modulations induced by peculiar motion in producing our signal, we cannot make quantitative statements on the luminosity distance bias. This is left for future work.
The trend in the FF shows that we start losing SNR for our GW190521-like binary at a distance ≈ 0.3 − 0.4 pc: this is where the orbital acceleration becomes non-negligible, in agreement with our previous results [27]. At these distances, a template including a -4PN term would be able to detect the source with higher SNR.
The trend in Fig. 7 indicates that we could not detect GW190521-like binaries around AGNs using pure vacuum GR templates for outer orbital radii below 0.3 pc. Binaries on tighter orbits could only be recovered using the modified model described in this manuscript, including the full Doppler and Shapiro modulations. The use of these templates would increase the size of the template-bank to scan in matchedfiltering by the number of independent outer orbit parameters (three), and thus increase the threshold required to claim a detection [48]. In practice, this means that we will most likely need to perform an "archival search" [82], that is, to re-analyze LISA data if a massive binary, not identified in the online analysis, is detected by the 3rd generation of groundbased GW detectors (such as Einstein Telescope and Cosmic Explorer). A ground-based detection would narrow down the search in the parameter space and allow for the use of modified (non-vacuum GR) templates, like the one presented here.

V. CONCLUSIONS
Interactions between BH binaries and the surrounding astrophysical environment are expected to have a negligible effect on the GW signals observable by LIGO and Virgo (although they may produce detectable electromagnetic counterparts [14]). The reason is that these environmental effects tend to become important at lower frequencies, where they have time to build up a more significant phase difference in the gravitational waveforms.
From this point of view, it is quite natural to focus our attention on stellar-origin BH binaries in the lower frequency (∼ mHz) LISA band. These systems are expected to be observable by LISA months or even years before they merge in the band of ground interferometers [83], and could probe their local environment, revealing the presence of gas in high densities or third bodies, possibly pointing at the formation mechanism of the binary itself or probing otherwise unexplorable galactic centers [20,27,31]. The same systems could provide important information on putative low-frequency, non-GR effects, such as BH "hair" and vacuum dipole emission [36,37]. It is therefore crucial to understand how matter and the environment affect these binaries, both because interesting in themselves and because they may bias tests of GR.
Recently, LIGO and Virgo have detected a particularly massive BH binary, GW190521, with total mass in excess of 100 M [3,4]. Such a large mass, as well as the possibility that this system may have an optical counterpart [14], may suggest that this binary formed in a gas-rich environment, e.g. an AGN disk, and that it might the first of many AGN (or GW190521-like) GW events.
In this work, we have extended our earlier analysis of GW190521-like binaries in LISA [27] to include two effects: the Doppler and Shapiro effects induced by the galactocentric motion of the binary. We proposed a way of producing accurate frequency-domain waveforms, combining the chirping and anti-chirping parts of the Doppler-modulated signal. By performing a full Bayesian analysis, we have found that the Doppler modulation and Shapiro time delay can induce significant effects on LISA waveforms if the binary is at distances of about a thousand gravitational radii (or less) from the AGN's supermassive BH. We also show that these effects permit breaking the degeneracy between purely vacuum wave-forms and ones including environmental effects, thus allowing LISA to measure the AGN parameters with high precision (galactocentric radius of the binary, central BH mass and AGN inclination angle). Not only does this have important implications for AGN astrophysics, but it will also allow for potentially disentangling environmental imprints on the waveforms from putative non-GR effects entering in the same lowfrequency band.
In our analysis, we focused on the accurate treatment of the leading effects of the third body on the waveform, by focusing on the induced time shifts. While subdominant, amplitude effects, neglected here, should eventually be included in the waveform model to avoid minor biases in luminosity distance measurements (see e.g. [23] and [84]) and to improve the determination of the environmental parameters.
Another natural extension of our work will be to consider galactocentric orbits highly aligned with the line of sight, for which the effect of gravitational lensing is non-negligible. Our accurate treatment of Doppler modulated signals will allow for a better assessment of the detectability of lensing in AGNs, which in turn might improve the precision with which LISA could measure the SMBH mass. Our method to treat the Doppler-induced signal modulation could also be used for systems on even tighter orbits around SMBHs, to better assess higher order effects relevant in this regime [31,39,73]. cussed in detail in Ref. [23] ) in the large separation limit. For simplicity, we neglect the cosmological background here. When the outer orbit is very wide, the standard SPA can be used to approximately compute the waveform in the frequency domain. Formally, this corresponds to expanding up to leading order in the separation of timescales between radiation reaction and outer orbital motion. We can then write [77,85], where we separated the vacuum signalh( f ) from the contribution due to the change in the time of arrival of the signal. At leading order the time-frequency correspondence is given by the SPA for the vacuum signal only 6 , with ϕ the Fourier domain phase of the vacuum signal. At leading post Newtonian (PN) order, where η is the symmetric mass ratio of the inner binary, we ignored a constant phase and a time shift, and v ≡ (πM f ) 1/3 . Notice that the frequency here is always the observed one, not that at the source.
Inserting Eq.(A3) in Eq.(A2) one gets At larger separations, we can expand the projected motion (3.8) in terms of the velocity and acceleration along the line of sight (pointing from the observer to the source) at t = 0, and ignore a constant r 0 , which defines the distance to the source in the vacuum waveformh. Then the observer-frame waveform can be written ash o where again we ignore terms that are constant and linear in the frequency, corresponding to a constant phase and a time shift. The first term corresponds to a Doppler shift and for a source at cosmological distance defines the redshift together with the usual cosmological contribution see also Ref. [23]. The second term is a -4PN correction due to the peculiar acceleration, and was explored in detail in Refs. [21][22][23].   The parameters are, from left to right: chirp mass, mass ratio, symmetric spin combination (χ PN ), asymmetric spin combination (χ − ), initial GW frequency, luminosity distance, azimuthal position of the observer in the source frame, ecliptic longitude and latitude, polarization phase, orbital angular velocity of the binary around AGN, projections of the orbital separation a • cos ι • , a • sin ι • , and initial centrogalactic orbital phase. The black lines mark the true parameters of the injected signal.