The peculiar acceleration of stellar-origin black hole binaries: measurement and biases with LISA

Nicola Tamanini, Antoine Klein, 3 Camille Bonvin, Enrico Barausse, 5, 6 and Chiara Caprini Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am Mühlenberg 1, 14476 Potsdam-Golm, Germany CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98 bis Bd Arago, 75014 Paris, France Institute for Gravitational Wave Astronomy and School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, United Kingdom Départment de Physique Théorique and Center for Astroparticle Physics, Université de Genève, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland SISSA, Via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste IFPU Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy Laboratoire Astroparticule et Cosmologie, CNRS UMR 7164, Université Paris-Diderot, 10 rue Alice Domon et Léonie Duquet, 75013 Paris, France (Dated: July 4, 2019)

We investigate the ability of the Laser Interferometer Space Antenna (LISA) to measure the center of mass acceleration of stellar-origin black hole binaries emitting gravitational waves. Our analysis is based on the idea that the acceleration of the center of mass induces a time variation in the redshift of the gravitational wave, which in turn modifies its waveform. We confirm that while the cosmological acceleration is too small to leave a detectable imprint on the gravitational waveforms observable by LISA, larger peculiar accelerations may be measurable for sufficiently long lived sources. We focus on stellar mass black hole binaries, which will be detectable at low frequencies by LISA and near coalescence by ground based detectors. These sources may have large peculiar accelerations, for instance, if they form in nuclear star clusters or in AGN accretion disks. If that is the case, we find that in an astrophysical population calibrated to the LIGO-Virgo observed merger rate, LISA will be able to measure the peculiar acceleration of a small but significant fraction of the events if the mission lifetime is extended beyond the nominal duration of 4 years. In this scenario LISA will be able to assess whether black hole binaries form close to galactic centers, particularly in AGN disks, and will thus help discriminate between different formation mechanisms. Although for a nominal 4 years LISA mission the peculiar acceleration effect cannot be measured, a consistent fraction of events may be biased by strong peculiar accelerations which, if present, may imprint large systematic errors on some waveform parameters. In particular, estimates of the luminosity distance could be strongly biased and consequently induce large systematic errors on LISA measurements of the Hubble constant with stellar mass black hole binaries.

I. INTRODUCTION
The long awaited detection [1][2][3][4] of gravitational waves (GWs) by the LIGO interferometer, followed by Virgo [5,6], opened a new era in the history of astronomy. Indeed, not only did these observations provide direct evidence for the existence of GWs (whose existence was previously assessed only indirectly by timing binary pulsar systems [7]) and for the validity of General Relativity (GR) in the highly relativistic strong field limit [8], but they also shed light on the physics of compact objects (black holes -BHs -and neutron stars -NS) [9,10]. Indeed, the coincident detection (to within 1.7 s) of the GW signal GW170817 and the gamma ray burst (GRB) 170817A has substantiated the long suspected connection between NS mergers and short GRBs [11]. Moreover, these joint observations have allowed constraints to be placed on the cornerstones of our understanding of gravity (Lorentz symmetry and the equivalence principle) [11], as well as on whole classes of gravitational theories modifying and/or extending GR (see e.g. [12,13]). GW observations have also started testing the very existence of BHs, down almost to the scale of the event horizon [14,15], where exotic non-GR physics is still possible, though tightly constrained [16][17][18]. Furthermore, as more and more BH detections accumulate, GW interferometers will reconstruct the mass function of these objects with increasing precision. This, together with future measurements of the eccentricity and spin of these binaries, will lead to better understanding of their astrophysical origin (see e.g. [19,20]).
Given this tumultuous succession of discoveries, it is hardly surprising that the European Space Agency (ESA) approved in June 2017 the Laser Interferometer Space Antenna (LISA) mission [21], whose launch is nominally scheduled for 2034 and which will target the mHz band of the GW spectrum. This band of frequencies is inaccessible from the ground, due to the seismic noise affecting the operations of terrestrial interferometers, and contains an impressive wealth of astrophysical GW sources. These include the merger of massive (∼ 10 4 -10 7 M ) BH binaries [22,23]; Galactic and extra-Galactic white dwarf binaries [24][25][26]; extreme mass ratio inspirals consisting of a stellar-mass BH and a massive BH [27,28]; and also the stellar-mass BH binaries that are targeted at higher frequencies by LIGO/Virgo [29]. Indeed, those binaries may be observable by LISA in the tens or even hundreds, especially if their mass is sufficiently large, for several years. These sources will then leave the LISA band before resurfacing in the LIGO/Virgo band, where they will merge [29].
While the signal-to-noise ratio (SNR) of these binaries in the LISA band will be at most of the order of a few tens or less (with most sources being indeed unresolvable individually) [29], these sources will accrue several inspiral cycles at the high end of the LISA sensitivity curve. This will potentially result in accurate determinations of the source parameters, such as BH masses and spins, eccentricity, sky position, etc [20,30]. This will also make these sources excellent probes of putative tiny deviations away from GR (orders of magnitude weaker than those detectable with ground interferometers) [31,32] or even environmental effects from the interaction of the binary with the surrounding matter [15,33].
Among these environmental effects, an especially attractive one is provided by the peculiar acceleration of the binary's center of mass (CoM) with respect to the observer. The peculiar velocity of the binary causes a shift in the signal's frequency, via the well-known Doppler effect. If the velocity is constant during the GW measurement, the Doppler frequency shift can be simply reabsorbed into a change of the (redshifted) chirp mass of the binary and into a shift of the luminosity distance. However, as shown in [34], if the velocity varies with time, i.e. if the binary's center of mass is accelerating with respect to the observer, the gravitational waveform changes in a (potentially) measurable fashion. This effect has a characteristic time (or frequency) dependence, and its amplitude is proportional to the magnitude of the acceleration component along the line of sight. The center of mass acceleration depends crucially on the binary's origin: formation in dense nuclear star clusters [35][36][37] or via accretion disks in active galactic nuclei (AGNs) [38][39][40][41][42] provides on average larger accelerations than, for instance, field formation [43].
In this paper, we investigate the conditions under which the acceleration can be measured by LISA, and its impact on the estimation of the binary parameters. The binary acceleration has two components: one of cosmological origin, due to the expansion of the Universe during the observation of the binary [44,45]; and the other of astrophysical origin, i.e. the peculiar acceleration caused by surrounding matter [34,46]. Our results confirm that the first component is not detectable with LISA [34], and we thus focus only on the second one 1 . We assume a constant peculiar acceleration on the timescale of the GW measurement, which is a good assumption for the astrophysical systems considered in this work. The magnitude and direction of the acceleration are kept as free parameters. By adopting the official LISA configuration of [21] and astrophysical binary BH populations calibrated to the LIGO/Virgo detection rate, we investigate via a Fisher matrix approach the precision with which the peculiar acceleration can be inferred from the data, improving on the estimates of [34,43]. We consider two values for the LISA mission duration: 4 years, the nominal duration, and 10 years, a very plausible extension [21]. We assess the improvement in the estimate of the acceleration due to the joint detection of the same event by ground based interferometers. Overall, our results confirm that LISA will be able to observe at least a few events with high enough peculiar acceleration, e.g. if the latter are produced in dense environments such as nuclear star clusters or AGN disk, and if the mission will last for 10 years. The detection of such peculiar acceleration by LISA would be extremely valuable, since it may provide information on the binary's formation mechanism, especially if the host galaxy can be identified electromagnetically with sufficient confidence (to identify for example the presence of an AGN or a nuclear star cluster). Note that the sky position of stellar mass BH binaries may be identified to within 10 deg 2 or better [30].
In addition, we also investigate the bias that the binary's acceleration would generate on the source parameters, if its effect is not properly accounted for in the GW waveform. Our analysis, which is based on a Fisher matrix approach and which is valid in the small-bias limit [47][48][49][50], extends and completes the investigation performed in [34], where the bias was calculated using simplified Monte Carlo simulations. We confirm that the bias is typically small if peculiar accelerations are below ∼ 10 −7 m/s 2 , but strong biases (larger than the inferred 1σ errors) may be present for binaries with larger peculiar accelerations, particularly those forming in AGN disks and nuclear star clusters, even for a 4 year LISA mission. Remarkably, the waveform parameter most strongly biased is the luminosity distance of the source. This suggests that strong peculiar accelerations can induce a large systematic error on measurements of the Hubble constant with stellar mass BHBs observed by LISA (see e.g. [30,51]).
The plan of the paper is the following. In the next section we derive the impact of the peculiar acceleration on the observed waveform associated with stellar mass BHBs. In section III we describe the procedure to implement the acceleration effect in the source catalogues and in the Fisher matrix code. In section IV we present our results and we conclude in section V.

II. THE EFFECT OF THE CENTRE OF MASS ACCELERATION ON THE GRAVITATIONAL WAVEFORM
In this section we show how the acceleration of the center of mass of a binary system emitting GWs affects the waveform signal detected by LISA. We follow the derivation presented in [34]. Note that, with respect to [34], we concentrate only on the redshift perturbation due to the homogeneous expansion of the universe and the peculiar velocity, which provide the dominant effects on the waveform. We therefore neglect the cosmological contributions due to the Bardeen potentials and the integrated Sachs Wolfe effect. Furthermore, we correct a sign error in the wave phase, c.f. Eq. (28), and add a new contribution in the wave amplitude, c.f. Eq. (33). In the following we provide the main points leading to the perturbed waveform in Fourier space, while the full derivation is reported in the Appendix.
We start by considering the GW signal in the source frame, which is given by [52,53] where M is the intrinsic total mass of the system, ν = m 1 m 2 /M 2 is the symmetric mass ratio, d C is the comoving distance from the source to the detector, n is a harmonic number, A n is the corresponding wave amplitude, φ src is the orbital phase (at the source), and y = (GM ω src /c 3 ) 1/3 is a post-Newtonian parameter used to describe the system, with ω src = 2πf src where f src is the orbital frequency of the system in the source frame. An observer at cosmological distances from the source would measure this signal redshifted by both the background expansion of the Universe and the presence of matter structures between the source and the observer, which modify the frequency at observer as follows where z is the redshift factor, which is a function of time.
The dominant contributions to the redshift are due to the cosmological expansion of the Universe and to the Doppler effect (see e.g. [54]) where a(t) is the scale factor, t src and t obs denote the local time of the source and of the observer, and v src = n · v src and v obs = n · v obs are their respective peculiar velocity along the line of sight n (here n points from the observer to the source). As shown in [34] the gravitational redshift and integrated Sachs-Wolfe contributions can be safely neglected. Moreover given that the phase does not change from the source to the observer, the time intervals are related as follows

A. Constant redshift
If z is assumed to be a constant, i.e. if both the expansion of the Universe and the peculiar velocity of the source are constant during the GW observation 2 , then the GW signal measured by the observer becomes simply A n e −inφ obs + c.c., (6) where M z = (1 + z)M is the redshifted total mass and d L = (1 + z)d C is the luminosity distance. The post-Newtonian parameter is unaffected because At leading post-Newtonian order, the frequency evolution equation in the source frame is given by and in the detector frame the same equation can be expressed as dω obs dt obs = 96νc 6 y 11 We can use this equation with the relation between the observed frequency ω obs and the PN parameter y to compute the time-frequency relation where t c denotes the coalescence time in the observer's frame and we have used that the frequency diverges at t c so that y(t c ) → ∞. Similarly, we can compute the observed orbital phase Note that the peculiar velocity of the observer may vary significantly during the time of observation, but this is automatically accounted for in GW analyses since the intrinsic motion of the detector is known with great precision.
where φ c = φ(t c ) denotes the phase at coalescence. We can then compute the Fourier domain version of this waveform using the stationary phase approximation (see e.g. [55][56][57][58]) h obs (f obs ) = (13) 2π n|φ obs | A n y 2 n e i(2πf obs tn−nφ obs −π/4) , where a dot denotes derivative with respect to time and all functions are evaluated at the stationary time t n for each n. We can first use Eq. (14) to find and then combine Eq. (13) with Eqs. (10)- (12) to obtaiñ For a constant redshift and at Newtonian order, this can be roughly regarded as the signal that is detected and analyzed by GW interferometers such as LISA and LIGO/Virgo.

B. Time dependent redshift
Since the expansion of the Universe and the peculiar velocity of the source are not exactly constant during the time of observation of the binary, the redshift acquires a small time-dependence. In this case, the waveform detected by an observer differs from the one derived above [34,44,45]. Let us briefly repeat the derivation presented in [34]. Under the assumption that both the cosmic expansion and the peculiar velocity are slowly varying over the observation time, we can linearly expand the redshift around a fixed chosen time t * obs . As shown in Appendix 1 (see also Eq. (36) in [34]), the redshift at t obs can then be written as where τ obs = t c − t obs is the observed time to coalescence and Y (z * ) is a small correction encoding the variation of the redshift between t obs and t * obs . For binaries that are observed close to the coalescence time t c , a convenient choice 3 is to take t * obs = t c . In this case, Eq. (18) reduces to where Here H(t) is the Hubble rate andv src andv obs are the line of sight peculiar acceleration of the source and the observer, respectively. Note that here and in the following we drop the subscript 'obs' and 'src' on the coalescence time t c , to simplify the notation. Quantities at the source are always evaluated at the coalescence time in the source rest frame t c src , whereas quantities at the observer are evaluated at the coalescence time in the observer rest frame t c obs . At zeroth order in Y c we can absorb the effect of the constant redshift z c by rewriting our expressions in terms of the redshifted mass M z = M (1 + z c ), leading again to Eqs. (16)- (17). In what follows we will perform all computations at linear order in Y c τ obs 1. We can now repeat the steps that led to Eqs. (16)-(17) for a time dependent redshift. At the lowest order Eq. (11) relates the observed time to coalescence τ obs to y. To compute the correction at first order, we can use the frequency evolution equation in the source frame, namely Eq. (9), to get the evolution of the PN parameter y, and then integrate to get the observed time τ obs = t c − t obs as where we have inserted Eq. (11) to perform the integral. The orbital phase is again given by Eq. (12), while the observed orbital frequency becomes To compute the effect on the phase of the harmonics of the Fourier domain waveform, we can simply write 2πf obs = nω obs (t n ) = n dφ obs dt obs (t n ) .
The time-frequency relation yields which inverted gives Then using this relation together with Eqs. (22) and (12), the phase (24) becomes This last equation (with n = 2) is also derived in Appendix 2 following a slightly different route. We can then introduce and we find We can also compute the relation from which we can infer an evolution equation for y f A time dependent redshift correction appears also in the amplitude of the GW signal. The details of this calculation are presented in Appendix 3. In terms of the frequency, the waveform amplitude in the restricted waveform approximation (i.e. the leading PN order amplitude of the second harmonic n = 2) is given bỹ where χ c ≡ c tc 0 dt/a(t) denotes the background comoving distance to the source evaluated at coalescence. The luminosity distance d L (z) is affected by perturbations, namely by gravitational lensing at high redshift and by peculiar velocities at low redshift (see eq. (A.31) and [59]). As shown in Appendix 3, these perturbations also affect the frequency-dependence of the amplitude. Indeed, the term in the square bracket of eq. (33) gets a contribution from the time variation of the luminosity distance, which is proportional to the time variation of the peculiar velocity. Note that the time variation of the gravitational lensing contribution is neglected here. This variation is proportional to the time variation of the gravitational potentials, which is expected to be much smaller than the time variation of the peculiar velocity 4 . Eq. (33) differs from the result in [34], which does not account for the time variation of the luminosity distance. For sources at low redshift, where χ c is small, the second term is the dominant one, whereas at high redshift the first one dominates.
The expressions for the phase (30) and the amplitude (33) have been derived under the assumption that the scale factor a(t) and the source peculiar velocity v src (t) do not vary too much between the window of observation and the coalescence time. This is a good approximation for the binaries we will consider in the rest of this work, since the physical mechanisms producing their peculiar accelerations act on time scales much longer than the observational period of LISA, in particular on galactic time-scales. Furthermore, as we will see in Sec. IV, the binaries for which the acceleration effect is most relevant are those that are observed close to coalescence. If other types of peculiar accelerations are analyzed, for example the perturbation of the binary's CoM due to a third close companion object, then the time variation of the acceleration must be taken into account and the analysis performed here will no longer be applicable (see e.g. [60][61][62][63][64]).

A. Implementation of the peculiar acceleration effect
Following the definition in [34], we parametrize the peculiar acceleration of the binary centre of mass with a parameter such thaṫ This form is motivated by the particular case of a binary describing a circular orbit around the galactic center. In this case, the centripetal acceleration takes indeed the simple formv where e is the direction of the acceleration and n the one of the line of sight, v src is the (circular) velocity and r is the distance from the galaxy center. The parameter corresponds to these quantities being normalized as follows: = v src 100 km/sec The argument above motivates the choice of pre-factor in Eq. (34); however, it is clear that can be used to parametrize any kind of peculiar acceleration, besides the centripetal one. Note that, in principle, the acceleration of the observer also plays a role (c.f. Eq. (20)). However, we assume that this is already accounted for and subtracted via the standard GW detection procedure. The parameter defined in (34) is convenient to parametrise the peculiar acceleration of the binary. We use it in the construction of the binary population catalogues, as presented in the following section. On the other hand, the effect on the gravitational waveform includes also the component due to the universe expansion (c.f. section II). When we implement the total effect of the centre of mass acceleration on the waveform, we do it in terms of a related parameter, α, defined in Eq. (39).
In section III C, we discuss the link among and α.

B. Astrophysical black hole binary populations
We simulate two different black hole binary (BHB) populations, each assuming a different individual mass distribution: the LogFlat model assuming a uniform distribution in log-mass dN/dm ∝ m −1 for each component of the binary, and the Salpeter model assuming a Salpeter mass function dN/dm ∝ m −2.35 , both with 5M < m i < 100M . LIGO/Virgo observations provide merger rate estimates in the local Universe for both of these models [3], that we assume to be R = 32 Gpc −3 yr −1 for the LogFlat model, and R = 103 Gpc −3 yr −1 for the Salpeter model. These median rates were recently revised to be R = 19 Gpc −3 yr −1 for the LogFlat model, and R = 57 Gpc −3 yr −1 for the Salpeter model [65], which would roughly halve the results presented in this paper. However, the median rates used here are close to the upper limits of the confidence intervals presented in [65], and thus can be considered optimistic results from this point of view.
For each of those models, we simulate the local black hole binary population as follows. We set a maximum comoving distance of d max C = 2 Gpc, a minimum initial orbital frequency of f min = 2 mHz, and a maximum initial orbital frequency of f max = 10 Hz (the initial frequency corresponds to the frequency at which the binary is emitting when LISA starts observing). Assuming a uniform distribution of the merger times, we can compute the corresponding initial frequency evolution assuming leading order PN frequency evolution: where C is a constant. We create six realisations of the binary black hole population, each for a different distribution of the acceleration parameter . We use a parameter E to denote the magnitude of the acceleration vector, and we create one population with = 0, and one family of distributions with E following a Gaussian distribution centered at E = E m with standard deviation E m . We use the following five values of E m = [10, 10 2 , 10 3 , 10 4 , 10 5 ] to create 5 distributions for E. These values have been chosen by equal logarithmic separation from 0 (no peculiar acceleration) to the highest expected accelerations for stellar mass BHBs living close to the galactic center ( = 10 5 corresponds to an acceleration of ∼ 10 −6 m/s 2 ). We then obtain the corresponding as the projection of the acceleration vector along the line of sight by multiplying E by a uniformly distributed number in [−1, 1]. In what follows we will denote these distributions as acceleration scenario 1 to 5, with scenario 0 being the distribution without any peculiar acceleration, i.e. where has been set to zero for each event.
We then compute the yearly merger rate inside a sphere of radius d max C centered on the Solar System based on the LIGO/Virgo estimates. For each realization, we then draw a Poisson distributed random number N P compatible with this expected rate. We then randomize the sky location and comoving distance assuming that binary systems are uniformly distributed inside the sphere, with an initial orbital frequency following a distribution proportional to f −11/3 , and individual masses distributed by the LogFlat or the Salpeter model. We randomize all additional vector components (spins, orbital orientation) uniformly distributed on the sphere, the initial orbital phase uniformly distributed inside [0, 2π], and the dimensionless spin magnitudes uniformly distributed inside [0, 1]. For each binary, we compute the merger time according to the time-frequency relation Eq. (22). We stop our simulation when the number of simulated systems with total mass M < 100M merging within one year, i.e. satisfying τ obs < 1 yr at the start of the mission, reaches N P .
For each value of E m and for both astrophysical mass distributions, we simulate 20 catalogues taking into account two different observational scenarios, corresponding to the LISA mission taking data for 4 or 10 years. We thus produce a grand total of 480 catalogues: 20 × 6 (values of ) × 2 (mass distribution) × 2 (LISA duration), spanning 3360 years worth of data. As we will see in section IV A, each catalogue contains a few hundred merging systems detectable by LISA, c.f. Table I.

C. Waveform generation and parameter estimation
For the present study, we use an inspiral-only Fourier-domain gravitational waveform including spinprecession, 2.5PN order harmonics, and 3.5PN phasing [66]. Note that all the quantities used in the following are in the observer frame, and thus we drop the "obs" subscripts from now on. We include the effects of the cosmic expansion and the peculiar acceleration by adding the linear effect in Eq. (32) into the frequency evolution equation in the following way: where we parametrize the acceleration effect by andω 0 stands for the standard 3.5PN frequency evolution equation without acceleration effects. Note that, neglecting the contribution of the homogeneous expansion of the Universe, which is always sub-dominant with respect to the peculiar acceleration for the systems analyzed in this work, one can relate α to by (40) As we will see, the relative errors on M z and ν measured by LISA are always much lower than the relative error on α, meaning that the latter dominates the relative error on . We therefore use the following approximation to estimate the measurement error on : This approximation holds especially in the relevant cases under analysis here, i.e. those for which the peculiar acceleration effect is sufficiently large to be measured by LISA.
Irrespective of the effect due to the modification to the phase evolution equation, the peculiar acceleration induces an additional effect on the GW amplitude, through a modification of the luminosity distance. It amounts to, as derived in Appendix 3, Since every waveform amplitude A n is proportional to 1/d L , using Eqs. (22), (31) and (39) we multiply our waveform by a factor 1 + 5α 128y 8 The Fourier transform of the waveform is then computed using a shifted uniform asymptotic (SUA) transform [66], where a dot denotes differentiation with respect to time and the coefficients a k,kmax are constants satisfying the following linear system of equations kmax k=1 a k,kmax + 1 2 a 0,kmax = 1, We simulated the LISA response with a low-frequency approximation, as presented in [67]. In order to evaluate the parameter estimation capabilities of LISA, we use a Fisher matrix analysis method (see e.g. [55][56][57][58]). We define the inner product in the space of signals as where a and b are waveform signals in the LISA detector, the tildes denote their Fourier transforms, the star denotes complex conjugation, and S n (f ) is the one-sided power spectral density of the instrument's noise. We use the LISA proposal noise curve [21,28] with a four-year galactic binary confusion noise: with acceleration noise, noise from other sources, and confusion noise For a given binary in our catalogues, we construct its waveform as described above as a function of its parameters h(θ), where θ is a 16-dimensional vector containing the 15 parameters describing a spinning compact object binary on a circular orbit, with the extra acceleration parameter α. The SNR ρ of this binary can be estimated as The Fisher information matrix for this binary is given by where ∂h/∂θ i is the derivative of the waveform with respect to the parameter θ i . We then construct its inverse, the covariance matrix Σ = Γ −1 . Its diagonal elements provide estimations of the statistical errors involved in the measurement of the parameters as We follow [47] to provide an estimate of the parameter estimation bias induced by neglecting the acceleration parameter α. We define the true waveform as where θ 15 is the parameter vector excluding the acceleration parameter α, and we define the approximate waveform as We then compute the 15-dimensional Fisher matrix Γ (15) using the approximate waveform h ap : and its inverse Σ 15 = Γ (15) −1 .
We then get an estimate of the bias that we would get in the estimation of the binary parameters by ignoring the acceleration effects as This method is valid in the high SNR limit (see [47] for a detailed discussion), but provides an efficient way of estimating the biases caused by a mismodelling of the intrinsic acceleration. Note that this method is accurate only if the biases induced by the mismodelling remain sufficiently small [50].

D. LISA mission simulations
For each of our 240 catalogues, we simulate two realizations of the LISA mission differing in the mission duration, that we take to be four years and ten years respectively. In each of these simulations, we compute the SNR for each binary. We separate the binaries into three distinct categories, depending on the time to merger τ c , defined as the interval of time between the beginning of LISA observation, denoted by t start obs , and the time of coalescence 5 : τ c ≡ t c − t start obs . Following [68], for τ c < 10 yr, we use an SNR threshold ρ thr = 15 for a LISA-only detection, and we assume an archival search from a groundbased detection with SNR threshold ρ thr = 9.5 for a multiband detection; for 10 yr < τ c < 100 yr, we assume a LISA-only detection with ρ thr = 15. For τ c > 100 yr, we assume a different search strategy is employed, leading to ρ thr = 10: Mangiagli et al. [69] found that a Newtonian phasing will be sufficient for unbiased parameter estimation if τ c > 100 yr. We use this as evidence that those systems will be morphologically similar to Galactic white dwarf binaries, for which searches have been implemented [70,71] and found an SNR threshold of ρ thr ∼ 10. Note that all of these studies neglected acceleration effects. However, in a realistic scenario, non-accelerating waveforms would be used for detection and accelerating waveforms would then be used for parameter estimation. Therefore, this justifies using these studies to infer an SNR detection threshold also in the accelerating case.
For the binaries with SNR above the detection threshold, we compute the expected statistical errors ∆θ i on the 16 parameters, as well as the expected measurement bias ∆ th θ i 15 on the 15 binary parameters excluding the acceleration effect. For the binaries that merge within 10 years since LISA starts taking data, we also simulate a coincident ground-based detection. We do so by removing the entry in the Fisher matrix corresponding to the time of coalescence parameter, which we assume to be measured exactly by ground-based observations, thus reducing the dimensionality of the Fisher matrix Γ. We therefore implicitly assume that binaries with ρ > 9.5 in LISA are certainly detected by whatever detectors are present on Earth if they merge in a time span of ten years since the beginning of the LISA mission.

IV. RESULTS
In this section we present the results of the analysis described above. We first study for how many events the acceleration effect is detectable and then we investigate the binary parameter space to identify regions where the effect is most likely observable. Finally we discuss whether the acceleration effect biases the measured value of the binary parameters.

A. Number of detections
The main question that we want to address here is the following: given the astrophysical BHB populations that we produced following the method described in Sec. III, for how many events the acceleration effect will be measurable? To answer this question we investigate the events for which the relative uncertainty on α, defined by ∆α/α where ∆α is obtained through Eq. (57), lies below 100%. Note that as mentioned after Eq. (40) the estimated relative error on α can roughly be regarded as the relative error on , especially when the peculiar acceleration is high so that ∆α/α < 1.
For the acceleration scenarios 0 to 3 we find no events for which ∆α/α < 1, irrespectively of the details of the astrophysical population (Salpeter or LogFlat) or of the LISA mission duration. This implies that the acceleration effect will not be detectable by LISA for moderate peculiar accelerations ( 10 3 ). This result is in agreement with previous studies [34,43] and it is here confirmed using realistic simulations of the astrophysical population of stellar-origin BHBs, over which we have imposed a distribution of moderate peculiar accelerations. Another immediate consequence of this result is the confirmation that LISA will not be able to measure directly the expansion of the Universe. This can be inferred because the effect is unobservable in scenario 0, where peculiar accelerations are set to zero and consequently only the cosmological acceleration is left, c.f. Eqs. (39) and (20). This is an expected result: Ref. [45] demonstrated that only more advanced space-borne GW detectors such as BBO or DECIGO might be able to directly observe the cosmic expansion through its effect on the gravitational waveform.
In the remaining two scenarios, namely 4 and 5, where the acceleration vector has a typical magnitude corresponding to = 10 4 and = 10 5 respectively, we find that the acceleration effect might be detected for some events. From Eq. (34) we see that these two scenarios correspond to a mean acceleration of respectively 10 −7 m/s 2 and 10 −6 m/s 2 , which can typically happen if the binaries form in dense astrophysical environments, for example in nuclear star clusters or AGN disks. In the particular case of a binary system orbiting circularly around a massive central BH with mass M BH at distance r, can be rewritten as In this case, = 10 4 would correspond for example to a binary system orbiting at 66 pc of a BH of mass M BH = 10 10 M (whose Schwarzschild radius is r s = 9.4 × 10 −4 pc).
In Table I we report the total number of LISA detections together with the number of events for which the estimated error on α is measured with a relative uncertainty below 100%, 50%, 30% and 10%. Results are shown for both a 4 years and 10 years LISA mission and for both BH populations: LogFlat (upper numbers in each line) and Salpeter (lower numbers in each line). We report detections by LISA only and for multiband events detected first by LISA and then by ground-based detectors such as LIGO and Virgo. For the latter class of detections, we select events with a time to coalescence (measured from the start of the LISA mission) below 10 years. The detection by ground-based interferometers is accounted for by fixing the value of the coalescence time t c to its true value, and it is practically implemented by eliminating the row and column corresponding to the parameter t c in the Fisher matrix.
From Table I it is clear that even for scenarios 4 and 5 the acceleration effect will not be measurable if we assume that the LISA mission will last 4 years. Only very rare events will have a peculiar acceleration strong enough to be barely detectable, implying that no definitive conclusions on the astrophysical processes of BHB formations will be extracted. The same conclusions can be drawn for the scenario 4 with a 10 years LISA mission, where the effect can be measured for few events only if a coincident detection with ground-based interferometers is made. Scenario 4 yields thus results similar to scenarios 0 to 3 where the number of events for which the acceleration effect is measurable is too low to extract any useful information on the BHB population.
The situation is different for scenario 5 if a 10 year LISA mission is considered. According to Table I LISA will be able to measure (with ∆α/α < 1) the acceleration effect in about 17% of the total detected events, and for more than about 4% of them (roughly ∼ 10 events) the peculiar acceleration will be determined with a relative precision better than 30%. These results drastically improve for multiband coincident detections: more than 40% of these events will present measurable peculiar accelerations, with about 16% and 3% of them presenting precisions better than 30% and 10%, respectively. Moreover, these percentages are independent of the BH population considered, since they roughly hold for both the Salpeter and LogFlat mass functions. Ta-LISA mission Acceleration LISA only LISA+Earth (tc < 10y) duration scenario Total 100% 50% 30% 10% Total 100% 50% 30% 10% TABLE I. Average number of events for which the parameter α is measured with a relative error below 100%, 50%, 30% and 10%. The upper numbers are for the LogFlat mass function, the lower ones for the Salpeter. "LISA only" denotes detections by LISA alone, while "LISA+Earth" are the detections for which the time of coalescence has been measured by ground-based interferometers, such as LIGO/Virgo (in this case an upper limit of 10 years has been applied on the time to coalescence from the beginning of the LISA mission). For LISA only, a detection threshold ρ thr = 15 is used if tc < 100 yr, and ρ thr = 10 otherwise; for LISA+Earth, a detection threshold of ρ thr = 9.5 is used.
ble I shows, however, that different choices for the BHB population model do influence the total number of events for which the acceleration is detectable, since the LogFlat mass function provides systematically higher event numbers, due to LISA being more sensitive to higher mass systems [29]. Note that the total number of detections quoted here are higher than estimated in some other works (e.g. [68]). This is due to our choice of a higher mass limit of 100M instead of 50M . These results show that a BHB population with high enough peculiar accelerations will imprint a distinguishable signature on the events detected by LISA only if the mission will reach its maximum possible duration. Such a population would be well representative of a characteristic BH formation channel where binaries are created in AGN disks close to the galactic center [43]. The results reported here thus support the science case for a LISA mission duration longer than the nominal 4 years period currently planned [21]. If LISA observes GWs somewhat continuously for almost 10 years, then we will have the chance to gather useful information on possible formation channels of BHBs as well as their distance from the center of the hosting galaxy (for BHBs orbiting a galaxy the peculiar acceleration is expected to be roughly proportional to the inverse of the distance from the galactic center, assuming simple circular motion).

B. Effect across the parameter space
The magnitude of the acceleration effect on the gravitational waveform depends on the BHB parameters: there are regions in the binary parameter space where the acceleration effect is stronger, and thus more likely to be detected. For example, from Eq. (38) one can roughly infer that, for the same value of α, the product M z ω should be maximized (c.f. Fig. 1). The question we want to address in this section is the following: for which values of the BHB parameters can we expect to measure the acceleration effect? In order to answer this question we focus on the only scenario where a sufficiently large number of detections is expected, which, according to Sec. IV A, corresponds to scenario 5 with a LISA mission of 10 years. However, irrespective of the LISA mission duration or the BH formation scenario, this analysis allows us to select, in a set of measured events, those for which it would be more plausible to detect the center of mass acceleration.
We focus our analysis on two parameters in particular: M z = M z ν 3/5 (the redshifted chirp mass), and f min (the GW frequency at the initial time of LISA observations). All other physically interesting parameters follow distributions roughly correlated with the SNR: the acceleration effect is better measured at lower distances, for symmetric mass ratios close to ν = 1/4 and for higher values of α. In Fig. 1 we show the probability of measuring the acceleration effect in the (f min , M z )-parameter space, where for each 2D bin this probability is defined as the ratio between the number of events for which the relative 1σ error on α is below 100% and the total number of detections in that bin. We report results for both mass functions (Salpeter and LogFlat) and with both LISA only detections and multiband detections involving ground-based observations. In each panel we also plot the lines to the right of which binaries will merge in less than 10 or 4 years after LISA starts taking data.
It is evident from Fig. 1 that the acceleration effect is measurable preferably for binaries with τ c 10 years. Note that multiband events are limited by the condition τ c < 10 years, so no events appear below the solid red line in the right panels of Fig. 1. We therefore find that optimal detection is obtained for times to coalescence close to the mission duration. This confirms what was obtained by previous studies (cf. Fig. 3 of [43]), which, however, did not consider simulations of realistic astrophysical BHB populations. This result can be understood by two balancing effects: on the one hand, the parameters of the binaries can be better measured when we are observing close to coalescence (very far from coalescence the waveform evolves indeed very slowly inducing strong degeneracies between the parameters of our waveform model). On the other hand, the dephasing generated by the binaries' acceleration accumulates over the time of observation, and therefore the longer we observe, the better α can be determined. The optimal situation to detect the acceleration effect is to observe almost up to coalescence.
We also plot on Fig. 1 the curve τ c = 4 years (dashedred), corresponding to binaries which coalesce 4 years after the beginning of observation. This curve roughly marks the upper edge of the parameter space region for which the acceleration becomes measurable for a consistent number of events with a LISA mission lasting 10 years. This explains why, with a 4 years mission, the acceleration effect is not detectable for many events: we need to observe a binary for more than 4 years before coalescence in order to reach a non-negligible detection probability. This plot also suggests that increasing the duration of the mission slightly above 4 years would already allow us to detect the acceleration effect with a 20-30 percent probability.
The results we obtain here, shown in Fig. 1, are coherent with what found in [43], despite the differences between the two analyses: in [43] there was no simulated binary population and the mass ratio and the redshift were fixed, furthermore both the waveforms and the LISA mission characteristics were different. The upper limit around redshifted chirp masses of 120 M and the lower limit around M z 10 − 20M are both defined by the lack of BHBs with higher and lower masses, respectively, in our population of events detected by LISA (the maximum chirp mass for our population is 87M in the source frame and ∼ 130M in the detector frame). It appears from our analysis that the peculiar acceleration is comparatively easier to measure for low mass binaries, in particular for LISA alone. We also observe that accounting for multiband observations substantially improves the detection probability. Note however that, due to the low number of events detected by LISA in a single mission lifetime in the region M z 20 − 40M , f min 12 mHz, more statistics (i.e. more simulated populations of BHBs) would be required in order to better characterize this region. Nevertheless, our results indicate that, when the real data is available, it will be worthwhile to check for the presence of peculiar acceleration in events having binary parameters close to the above mentioned region, as they are more likely to be affected by the acceleration effect.

C. Biases
We finally investigate the bias induced by the acceleration effect in case it is neglected in the parameter estimation analysis. The questions we want to address here are the following: how many events will be biased by the acceleration effect if its contribution to the GW waveform is ignored? Which waveform parameters will be biased the most and for which values? To address these questions we estimate the bias on each waveform parameter of each event using the bias estimator discussed in Sec. III, namely Eq. (61). We analyze the results derived with LISA in this section and we do not consider coincident detections with Earth-based interferometers. We are in fact interested in understanding how LISA parameter estimation might be biased by the peculiar acceleration effect and how this might affect the expected coincident detection with ground-based instruments (e.g. with regards to the expected time of coalescence or sky position).
In Table II we report the percentage of biased events for the scenarios with non-zero peculiar acceleration. We average over the 20 realizations we have constructed, considering a LISA mission durations of 4 and 10 years, the acceleration scenarios 1 to 5, and both BH mass functions: LogFlat (upper numbers) and Salpeter (lower numbers in each line). For each waveform parameter an  Table II is defined as the number of events for which at least one parameter is biased over the total number of events.
From Table II we first notice that for each parameter the fraction of biased events is roughly the same for both BH populations (the Salpeter mass function systematically provides slightly lower percentages). This implies that the bias results of our analysis are roughly independent of the details of the underlying BHB astrophysical population. We recall though that the number of detections for the LogFlat mass function is higher than the Salpeter case (cf. Table I), meaning that in the latter case we would expect a lower absolute number of biased events.
From the last row of Table II it is furthermore clear that few if any events will be biased in acceleration scenario 1, irrespectively of the BH mass function (Salpeter or LogFlat) and of the LISA mission duration. The situation changes for acceleration scenarios with higher peculiar accelerations. For the acceleration scenario 2, we find that ∼ 14% of events will be biased if the LISA mission lasts 10 years, and a negligible amount if it lasts 4 years. For acceleration scenario 3 we find that ∼ 40% of the events are biased if the LISA mission lasts 10 years, and about 16% if it lasts 4 years. For acceleration scenarios 4 and 5, most events will be biased if the mission lasts 10 years, up to ∼ 96% for scenario 5.
From Table I we learned that a non-negligible fraction of events of scenario 5 with a 10-year LISA presents measurable peculiar acceleration. Therefore, finding many biased events in this case is expected. The acceleration effect will be measurable for some of the biased events, meaning that, for these, the bias on the waveform parameters may be corrected using the right, accelerated waveform. Somehow more surprising is the large number of biased events obtained in scenarios with smaller acceleration. According to Table I, only very few events have measurable peculiar acceleration in these cases, while here we find that a significant fraction of events (reaching ∼ 75% for acceleration scenario 5 and a 4-year mission, and scenario 4 and a 10-year mission), at least one of the waveform parameters is biased. One possible explanation is that the parameter α is degenerate with the other parameters. In this case, α is not large enough to break this degeneracy and allow for a detection of the effect. However, when α is not included in the modelling of the signal, the effect is reabsorbed into a shift in the other parameters, particularly the ones that are degenerated with α. As a consequence, in this scenario the acceleration effect is not detectable, but the bias induced on the other parameters can be important.
It is particularly important to remark that the parameter estimation of BHBs formed in AGN disks (acceleration scenario 5) might be biased by the peculiar acceleration effect even for the LISA nominal mission duration of 4 years. Note that in this case, adding α to the modelling of the signal would remove the bias on the other parameters. This is the case even for the events for which α is not measurable, at the expense of increasing the other  parameters errors.
We can now look at which waveform parameters are more biased on average. From Table II it appears that the parameter which is biased in the highest percentage of cases is the luminosity distance d L . We point out that this result is particularly important for LISA applications of stellar mass BHBs as standard sirens [30,51], since a possible systematic error on d L will negatively affect the accuracy with which cosmological parameters can be measured. In fact since at low redshift a bias of 1σ on d L directly translates into a bias of 1σ on the Hubble parameter H 0 . In Table III, we show the percentage of events for which the bias on the luminosity distance is larger than the 1, 2 and 3σ error on d L , respectively. We see that more than ∼30% of events will present a bias bigger than 3σ in the acceleration scenario 5 for 4 years of LISA observations. Similarly, in the acceleration scenario 4 for 10 years of LISA observations we find that this number is ∼25%. This implies that in these scenarios any measurement of H 0 can be strongly affected by the acceleration effect, introducing a systematic error in low-redshift cosmological measurement with LISA. The situation is even worse in the acceleration scenario 5 for 10 years of LISA measurement, where up to ∼80% of events will present a bias in d L higher than its 3σ error. The solution to avoid this bias is then to introduce the acceleration effect in the modelling of the waveform. Even if the acceleration parameter α is not large enough to be detected, by adding its effect to the waveform we ensure that the other parameters, including the luminosity distance, are not biased. This of course comes at the cost of widening the error bars, especially if α is degenerated with other parameters. Fortunately the sky localization of the source, represented by the angles θ N and φ N , seems not to be substantially biased, except in the most extreme case (10-year mission with acceleration scenario 5), where we expect to detect the acceleration effect for a large number of events (cf. Table I).
This is good news not only for standard sirens 6 , but also for multiband and multi-messenger detections (although counterparts are not really expected for BHB). Similarly to what happens for the luminosity distance, the direction of the orbital angular momentum of the binary, represented by the angles θ L and φ L , is also biased for a large fraction of events in the highest acceleration scenarios. On the other hand, the bias on the time of coalescence t c and the phase at coalescence φ c appears at substantially lower values of the peculiar acceleration. Small fractions of events are biased even for acceleration scenarios 2 and 3. These two parameters seem therefore to be particularly sensitive to the acceleration effect, a result to which [34] already hinted 7 . A possible interpretation of this fact is that t c and φ c are more degenerate with α than the other binary parameters.
We can now turn our attention to the parameter space, and investigate for which values of the binary's parameter a high fraction of detected events is biased. We again focus only on the parameters M z (redshifted chirp mass) and f min (initial frequency), to better compare with the analysis of Sec. IV B. In Fig. 2, we plot the probability of being biased, defined in each 2D bin as the ratio between the total number of biased events (those with bias larger than the 1σ error, as done for Table II) and the total events in the bin, in the (M z -f min ) parameter space. We only show the result for acceleration scenario 5, a 4-year LISA, and both Salpeter and LogFlat mass functions: this represents the most interesting case for a LISA mission of nominal duration. Indeed, the same figure in the case of a 10-year LISA mnission is less informative, as the results are biased across almost all the parameter space (cf. the numbers in Table II). Fig. 2 shows that events with the time to coalescence close to and lower than the mission duration have higher probability of being biased. In agreement with the results of Sec. IV B, for these events we expect relevant systematic effects induced by the peculiar acceleration. Moreover, as already hinted by Table II, many events for which the acceleration effect is not measurable turn out to be biased: e.g. events with τ c larger than the mission lifetime, as shown in Fig. 2. Of particular interest for multiband detections are low mass events, as well as events with τ c shorter than the mission duration (events lying above the red curve in Fig. 2), as they merge at higher frequencies and in a relatively short time. These events appear to have a higher probability of being biased. 6 A good sky localization is required to apply the so-called "statistical" method for standard sirens, where the sky position of GW sources is cross-correlated with galaxy catalogues [72][73][74]. 7 Note, however, that the analysis of [34] concentrated on very low-mass binaries, and on the probability of recovering tc in less than one minute. Since low-mass binaries are extremely far from coalescence, big absolute errors on tc are not at all surprising.

D. Intermediate mass black holes
We have also checked the detectability of the peculiar acceleration for intermediate mass BHBs (IMBHBs). By generalizing the procedure described above, we have constructed a population of IMBHBs with single BH masses up to 500 M . Since the rates of IMBHBs are presently unknown and poorly constrained, we have not used the normalization of such a population for our analysis, and we have abstained from making any prediction regarding the number of detections that LISA could attain. The main finding of this analysis is that, in analogy with stellar-origin BHBs, IMBHBs whose time to coalescence at the beginning of observations matches the mission duration, present a higher probability of yielding a measurable peculiar acceleration effect. Although the merger rate of an intermediate mass population of BHBs is completely unknown, it is interesting to note that they might actually be efficiently created in AGN disks [38,75,76], namely our scenario 5. The astrophysical models leading to IMBHs in AGN disks may thus be efficiently tested by LISA, which will certainly detect these systems and measure their peculiar acceleration if the mission lifetime will be extended over its nominal 4 year duration.

V. DISCUSSION AND CONCLUSION
The present work has been dedicated to assess the potential of LISA to detect the peculiar acceleration of stellar-origin BHBs. Based on the idea that the peculiar acceleration of a binary alters its gravitational waveform [34], we investigated how well it can be measured, using simulated BHB populations.
Our results confirm that the BHBs which are more likely to present a measurable peculiar acceleration are those that will merge close to the end of the mission [43]. The reason for this is twofold: on the one hand, the dephasing due to the peculiar acceleration accumulates with the time of observation; on the other hand, the uncertainty on the other waveform parameters decreases when observed close to coalescence, as many degeneracies in the phase are broken by the fast frequency evolution of the waveform.
Furthermore, we find that with the nominal LISA mission lifetime (4 years), no significant measurement is possible, irrespective of the magnitude of the peculiar acceleration considered (within plausible scenarios). On the other hand, we show that a LISA mission lasting 10 years will be able to detect the peculiar acceleration of a relevant fraction of binaries, if these form in dense astrophysical environments, for example in nuclear star clusters or AGN disks, where peculiar accelerations might reach values of ∼10 −6 m/s 2 (corresponding to ∼ 10 5 ; cf. Eq. (34)). Our results also suggest that, increasing the mission duration to little over 4 years might already be enough to discriminate between different BHB formation channels, shedding new light on the processes at work when BHBs are created. The peculiar acceleration effect investigated here will moreover be complementary to other effects expected from stellar BHBs in close orbit around a massive BH, in particular large eccentricities up to merger due to Eccentric Kozai-Lidov evolution [77]. A LISA joint detection of systems with observable peculiar accelerations and large eccentricities would constitute a strong evidence for BHBs forming in AGN disks or globular clusters.
Low mass binaries with a relatively short time to coalescence (we put a cutoff of 10 years from the start of LISA observation) can also be detected by Earth-based GW detectors at merger. We therefore also investigated how the acceleration effect may be measured in the case of joint LISA/Earth-based detection. We find that the number of events for which the effect becomes observable increases, reaching about 40% of all detected events in the highest acceleration case. As expected, the average uncertainty on the peculiar acceleration also improves for these joint observations, suggesting that the measurement of peculiar accelerations may constitute an interesting science case for multiband searches between LISA and Earth-based interferometers such as LIGO and Virgo.
For our analyses we assumed a LISA noise sensitivity curve as given in [21]. The high frequency part of this curve is based on a single link optical measurement system (OMS) noise of 10 pm/ √ Hz. Note that in the LISA Science Requirements Document [78] a margin has been inserted on this noise contribution, increasing it up to 15 pm/ √ Hz. Since our results are based on LISA observations at the high end frequencies of its sensitivity curve, the possible loss of sensitivity induced by this margin would affect them. As a rough estimate we are expecting the sensitivity to worsen by a factor of 1.5, and thus SNRs and measurement errors to degrade by the same factor. This would entail the loss of about 30% (50%) of the stellar-mass BHBs detections, and about 55% (22%) of the events with acceleration measurements with error less than 100%, for LISA-only (coincident) detections. Nevertheless, the number of events with acceleration measured to better than 10% in coincident LISAground based detections remains unchanged (such events are characterised by high SNR, thus the loss in sensitivity is less severe on their measurement).
We also studied whether the waveform parameters would be biased if the acceleration effect is ignored in the data analysis procedure. We computed a simple estimate of the bias on each waveform parameter, and studied how many events would experience a non-negligible bias. Our results clearly show that binaries possessing rather large peculiar accelerations (scenarios 3 to 5) will be biased even for a LISA mission duration of 4 years. We also found that binaries with a time to coalescence similar to the LISA mission lifetime will have a higher probability of being biased. The bias therefore follows the same parameter dependence of the measurement of the acceleration effect, but many events for which the peculiar acceleration would not be measurable are instead biased. This suggests caution when performing data analysis in the parameter space of stellar-mass BHB systems coalescing close to the LISA lifetime, since for these systems, ignoring the peculiar acceleration could lead to wrong estimates of the BHB physical parameters.
Furthermore, we have found that the parameter most often biased is the luminosity distance of the system, which is a fundamental quantity for cosmological applications of GW observations, in particular standard siren analyses. If the peculiar acceleration is not correctly taken into account, large systematics could affect the estimates of the cosmological parameters, e.g. H 0 , derived from LISA stellar-origin BHBs. The time to coalescence can also be biased: this occurs for a smaller fraction of events than the luminosity distance, but for lower values of the peculiar acceleration, suggesting caution also when estimating the statistical error on this parameter without inserting the acceleration in the waveform. If present, such a bias might affect multi-messenger and multiband searches which rely on an accurate estimation of the merger time of BHBs.
Finally, it is also interesting to note that the peculiar acceleration effect induces a perturbation on the phase of the gravitational waveform which has the same -4PN dependence in frequency of the effect given by a possible running of the gravitational constant [79]. In particular, if we compare the effect on the waveform phase given bẏ G (namely the time variation of the gravitational constant G) as computed in [79], with the acceleration parameter α given in Eq. (39), we finḋ This result implies that studies of the running of the gravitational constant with GWs, and in particular with stellar mass BHBs, might be strongly affected by the possible presence of large peculiar accelerations. The investigation of these issues is left for future works.

ACKNOWLEDGMENTS
This project has received funding (to E. B.) from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. GRAMS-815673; project title "GRavity from Astrophysical to Microscopic Scales"). This work was also supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. C. B. acknowledges support by the Swiss National Science Foundation. This work was supported by the center National d'Études Spatiales, by the COST Action CA15117 "Cosmology and Astrophysics Network for Theoretical Advances and Training Actions" through a Short Term Scientific Mission (STSM), and by the COST Action CA16104 "Gravitational waves, black holes and fundamental physics" (GWverse). COST actions are supported by COST (European Cooperation in Science and Technology).

Appendix: Derivation of the acceleration effect in the gravitational waveform
In this appendix we show in detail how the gravitational waveform produced by a binary inspiral is perturbed by the acceleration of the CoM of the binary system. We reproduce the results of [34] but correct a sign mistake in the wave phase and add a contribution in the wave amplitude.

Derivation of the redshift perturbation
We first derive the time dependence of the redshift for a GW source, expanding over suitable small quantities. We start from the general definition of the redshift as measured by a cosmological observer in an FRW universe (see e.g. [54]) where v src = n · v src and v obs = n · v obs are the peculiar velocities of the source and the observer (i.e. the velocity that the objects have with respect to the Hubble flow), respectively, along the line of sight (identified with the unit vector n), a(t) is the scale factor and we explicitly write all dependences on the source local time t src and the observer local time t obs . Here we ignore the contributions from gravitational redshift and integrated Sachs-Wolfe since, as shown in [34], they are sub-dominant. We consider v src /c 1 and v obs /c 1 as the non-relativistic astrophysical peculiar velocities of the source and the observer. Following [34] we want to expand the time dependent quantities in the r.h.s. of Eq. (A.1), which are assumed to be slowly varying in time, around a fixed chosen time in the source rest frame t * src and its corresponding time in the observer rest frame t * obs . In other words t obs = t * obs + (t obs − t * obs ) = t * obs + δt obs , where the time interval δt obs is such that the functions in Eq. (A.1) vary slowly during this interval. Similarly, we can expand t src around t * src , using δt src = t src − t * src as a small time interval. In what follows, we will expand the timedependent functions linearly in δt obs and δt src , and keep only the terms at first order in v src /c and v obs /c. Given the assumptions above, Eq. (A.1) can be expanded as 1 + z a(t * obs ) +ȧ(t * obs )(t obs − t * obs ) a(t * src ) +ȧ(t * src )(t src − t * src ) At zeroth order in the time evolution this expression reads 1 + z a(t * obs ) a(t * src ) 3) where we define z * = z(t * obs , t * src ) as the redshift at the time t * . Note that z * contains both the background expansion and the Doppler contribution. We can use this relation at the lowest order to express (t src −t * src ) in terms of (t obs − t * obs ): where we have defined −v obs (t * obs ) c . (A.6) We want to rewrite t obs − t * obs in terms of the time to coalescence τ obs = t c − t obs , where t c is the time at coalescence. By noticing that t obs − t * obs = t obs − t c + t c − t * obs = −τ obs + τ * obs , (A.7) we obtain 1 + z = (1 + z * ) 1 − 2Y (z * )(τ obs − τ * obs ) . (A. 8) For binaries that are observed close to the coalescence time t c , a convenient choice is to take t * obs = t c which yields τ * obs = 0 and consequently Eq. (19). Note that another possibility would be to choose t * obs as the time when LISA starts taking data t * obs = t st . For the binary systems analyzed in this work, τ st = t c − t st is at most of the order of thousands of years, so that τ st Y (z st ) 1. At first order we can therefore rewrite eq. (A.8) as Comparing this with eq. (19) we see that the expressions are equivalent at lowest order in Y τ obs . Indeed the constant redshift z c andz st are not observable, since they are degenerated with the mass, and therefore the difference between them does not matter. Furthermore the difference between Y c and Y st is of second order in the expansion, since it is due to the evolution of Y between the time when we start observing and the coalescence time. Therefore, for the binaries we are interested in, one can equivalently set t * obs = t st or t * obs = t c .

Derivation of the perturbed phase
Now that we know how the redshift changes over the time of observation of the binary, we want to understand how this affects the gravitational waveform as measured by the observer. We start by deriving first the effect in the phase and then we will consider how the amplitude becomes perturbed.
The GW frequency evolution equation at the source is given by [52,53] df src (t src ) dt src = 96 5 π 8/3 GM c Integrating this expression between t obs and t c (time of coalescence) and using that the frequency diverges at t c , so that g(t c ) → ∞, one obtains To this, we need to add the contribution from the luminosity distance in Eq. (A. 25). As shown in [59], the luminosity distance is affected by matter perturbations. At high redshift the dominant contribution is the one from gravitational lensing, whereas at low redshift it is due to the source peculiar velocity d L (z) =(1 + z(t src ))χ(t src ) 1+ (A.31) where χ(t src ) = c tsrc 0 dt/a denotes the background comoving distance to the source, ∆ Ω is the transverse Laplacian and φ and ψ are the two metric potentials. As before, we expand quantities at t src around the time at coalescence t c . We neglect in the following the contributions coming from the background evolution, since they are strongly subdominant with respect to the source peculiar acceleration. We neglect also the contribution due to the evolution of gravitational lensing. This contribution would indeed be proportional to the time variation of the metric potentials,ψ andφ. In a matter dominated universe, the potentials are constant in time, and these terms are exactly zero. The presence of a cosmological constant does induce a time variation of the potentials, but this time variation is much smaller than the time variation of the peculiar velocity and can therefore be safely neglected. We obtain then d L (z) = (1 + z c )(1 − 2Y c τ obs )χ c (A.32) where in the last line we have neglected terms that are second order in the perturbations. Including this expression into (A.25) we obtain A c (f obs ) = 4 d L (z c )c 4 (GM z ) 5/3 (πf obs ) 2/3 (A.33)  which agrees with Eq. (33) reported in the main body of the paper (recall that M z = ν 3/5 M z ). Note that the acceleration correction in the amplitude computed here differs from the one reported in [34] which didn't account for the velocity contribution in the luminosity distance. In all the analyses of this paper we took into account this amplitude perturbation by computing H c χ c assuming a standard ΛCDM cosmology.