Testing gravitational wave propagation with multiband detections

Effective field theories (EFT) of dark energy (DE) -- built to parameterise the properties of DE in an agnostic manner -- are severely constrained by measurements of the propagation speed of gravitational waves (GW). However, GW frequencies probed by ground-based interferometers lie around the typical strong coupling scale of the EFT, and it is likely that the effective description breaks down before even reaching that scale. We discuss how this leaves the possibility that an appropriate ultraviolet completion of DE scenarios, valid at scales beyond an EFT description, can avoid present constraints on the GW speed. Instead, additional constraints in the lower frequency LISA band would be harder to escape, since the energies involved are orders of magnitude lower. By implementing a method based on GW multiband detections, we show indeed that a single joint observation of a GW150914-like event by LISA and a terrestrial interferometer would allow one to constrain the speed of light and gravitons to match to within $10^{-15}$. Multiband GW observations can therefore firmly constrain scenarios based on the EFT of DE, in a robust and unambiguous way.


I. INTRODUCTION
Identifying the source of the current accelerated expansion of the Universe represents one of the biggest challenges in cosmology today. Models ascribing cosmological observations to a cosmological constant or to a dark energy (DE) component have been proposed, but are hardly natural from a microscopic point of view [1]. Given that all evidence for cosmological acceleration is gravitational in nature, it is possible that cosmological data can be explained without any such components, if deviations from General Relativity (GR) occur in the low energy regime, relevant for cosmology.
The simplest option for modifying GR is to include an extra scalar degree of freedom, non-minimally coupled with gravity, and whose time-dependent profile causes cosmic acceleration with no need of a cosmological constant. Over the past decades several scalar-tensor theories have been proposed, with increasing degrees of generality, from Fierz-Jordan-Brans-Dicke [2] theory, to Horndeski [3], beyond Horndeski [4,5], to DHOST [6][7][8]. A powerful theoretical framework, the effective field theory (EFT) of DE (see e.g. [9][10][11][12]) allows one to encapsulate the predictions of such scalar-tensor theories in terms of a set of operators, describing the possible couplings of the scalar field with gravity. Such an approach is useful for investigating the regime of validity of a theory and its energy cutoff, and for confronting scalar-tensor theories with cosmological data in terms of a small number of parameters: see [13] for a comprehensive review.
The operators appearing in the EFT of DE generically include non-minimal, derivative couplings of the metric to the scalar field. Once the scalar acquires a timedependent profile -as expected if it drives cosmic selfacceleration -the resulting speed of gravitational waves (GW) is different from the speed of light.
Using this fact, significant constraints on scalar-tensor theories have been placed by the detection of GWs in the LIGO-Virgo-KAGRA (LVK) band. In particular, the measured times of arrival of the GW signal from the neutron star (NS) merger GW170817, and of the electromagnetic signal from GRB 170817A, were coincident within two seconds [14][15][16]. This observation allowed the relative speeds of GW and gamma rays to be constrained to within a fractional difference 10 −15 [14], which in turn sets tight bounds on existing gravity theories [17][18][19][20][21][22][23][24][25]. However, the frequency of the GW170817 signal is expected to be above the cutoff in these EFTs of DE. Hence, it may be natural for the speed of GWs to be luminal at LVK frequencies while still departing from luminality at much lower frequencies [26].
Thanks to its lower frequencies of operation (in the mHz band), the future space-borne Laser Interferometer Space Antenna (LISA) [27] will provide GW observations at lower energies than the EFT strong coupling scale, thus potentially constraining scalar-tensor theories in a more robust and unambiguous way. However, the existence of electromagnetic counterparts for LISA is highly uncertain [28,29]; hence bounding the propagation speed of GWs with the massive black hole equivalent of standard sirens may not be feasible in practice.
As we discuss in this work, a promising alternative may instead be provided by multiband observations of stellarorigin black hole binaries (SOBHBs) with LISA and LVK (or a third generation terrestial GW detector). In fact, Ref. [30] shows that GW signals from this family of binaries, already detected by LVK, will be also observable with LISA in their early-stage, low frequency inspiral. LISA observations will take place months to years prior to the final merger in the band of terrestrial detectors, allowing for the merger time to be predicted to within tens of seconds 1 . For this reason, Refs. [31,35] proposed to use multiband observations of SOBHBs as a test for the presence of non-GR gravitons, which would lead to additional GW emission from the binary (on top of the usual tensor emission described at leading order by the quadrupole formula), and thus to an earlier merger. Additionally, [36][37][38] showed that multiband GW observations can improve theory-agnostic bounds on deviations from GR in the inspiral as well as test specific alternative theories of gravity, while [39,40] emphasized how multiband detections allow for breaking degeneracies and performing multiparameter tests of theories beyond GR.
In this paper, we show that a SOBHB merger event, which occurs in the band of terrestrial detectors at a time different than that predicted from LISA observations, may also be caused by a non-linear graviton dispersion relation. We assume a scenario whereby GW generation is described by GR (possibly thanks to non-linear "screening mechanisms" [41][42][43][44][45][46], see e.g. [47,48] for reviews) and only GW propagation is affected by modified gravity. In particular, we make the hypothesis that the propagation speed of GW is c T = c EM (where the EM stands for Electro Magnetic) in the band of ground-based experiments (as a result of new physics near the strong coupling scale of the EFT), while c LISA = c EM − δc in the LISA band, with δc/c EM very small. Then the LISA GW emission would be described by the same strain signal h as in GR, but evaluated at the retarded time t − D/c LISA (where D is the distance of the source; notice that since SOBHBs are only observable at low z 0.1, cosmological effects on time and distance can be safely neglected). A LISA observer would then predict the merger to occur at a time t LISA = t p + D/c LISA ± σ t , where t p is 1 These results have been confirmed by full Bayesian analyses [31,32], at least for quasi-circular binaries. More recently, Ref. [33] found much larger uncertainties on the predicted merger time than Refs. [31,32]. However, the disagreement disappears [34] for quasi-circular binaries, once a higher-PN order (and thus more accurate) expression for the merger time is employed in the analysis of Ref. [33].
the time at which the strain amplitude h(t) is the largest (which corresponds to coalescence), while σ t ∼ 5 sec is the statistical error on the merger time estimate coming from the LISA data analysis. However, when later detected in the LVK band, the merger actually happens at t Earth = t p + D/c EM , i.e. at a different time than predicted by LISA. Since the merger time in ground-based detectors can be measured with a precision of milliseconds, it is possible to constrain where, as mentioned above, we assume δc/c EM very small. For the distance of the first observed GW event, GW150914, D ≈ 410 Mpc and this bound becomes |δc|/c EM 10 −16 at 1σ confidence level (see Appendix A for more details). In this paper we explore this idea, applying a Fisher Information Matrix (FIM) forecast to LISA data analysis, and showing that the estimate of Eq. (1) holds. Multiband observations of SOBHBs with LISA can therefore severely constrain (or confirm!) scenarios based on the EFT of DE, in a more robust and unambiguous manner than ground-based observations alone. First, though, we start with explaining our theoretical motivations for a frequency-dependent speed of GW. During the last parts of this analysis, we became aware of similar methods and results being derived in [49]. The precise details of that analysis differ from the current work in places: for example, the authors of [49] consider a 'distinguishability criterion' to assess if the time-shift of their waveforms induced by varying c T is recognisable. In the present work we instead perform a FIM analysis, which allows us to make statements about simultaneous constraints on regular and modified gravity parameters, and their degeneracies. The conclusions presented here are nevertheless in perfect agreement with those of [49].
We work in units where = 1 and where the absolute speed of massless and free fields in the vacuum is c = 1. Note however that we keep the speed of light denoted by c EM , as in principle this could be slightly different. The Planck scale M Pl is related to Newton's constant by M Pl = (8πG) −1/2 .

II. THEORETICAL MOTIVATIONS
For sake of concreteness, we present an example of a low-energy EFT, where the speed c T of GWs is a function of frequency. We emphasize however that the methods and conclusions derived in this work apply well beyond the scope of this illustrative example. As a warm-up, we first consider the EFT for a single scalar field φ living in flat spacetime with low-energy interactions of the form where G a function depending on the precise details of the high-energy (UV) completion, and Λ the cutoff. This low-energy EFT may be obtained by integrating out massive modes H of mass above Λ, with interactions involving H(∂φ) 2 as described in [26]. We are not interested in the precise details of the UV completion here, but we expect the UV limit to be insensitive to the spontaneous breaking of Lorentz invariance induced at low energies by a scalar background. In fact, as mentioned above, if φ were to play the role of DE field causing self-acceleration, it would acquire a time-like profile. We can parameterize it as where t the physical time, and α a dimensionless constant that may be used to relate the scale of the background with that of the low-energy EFT. Fluctuations on top of the background scalar in eq. (3) are characterized by a frequency-dependent phase velocity c s (k), given by In [26] a power-law profile for G was considered. However, in principle the function G can potentially be sharper, and describe an abrupt transition between a constant non-luminal low-energy speed at low-energy c s (k Λ) ∼ 1 − δc, and a luminal speed at higher frequencies c s (k Λ) = 1.
Turning our attention back to GWs, derivative mixing terms between gravity and other light degrees of freedom are motivated by the EFT of DE, as explained in Section I. For instance, in the presence of a light scalar field φ, the Einstein-Hilbert action may be supplemented with non-minimal mixings while the EM sector remains minimally coupled to gravity, where generic dark-energy EFT operators are included within the quantity G µν given by In principle one could also expect the EFT to include pure-potential terms but those are irrelevant to the current analysis and would be absent or suppressed if the scalar is endowed by a shift symmetry. In the previous equation, we can in principle consider any covariant contraction of the Riemann tensor R, with an arbitrary number of higher derivatives. It is convenient, though, to start the expansion with Horndeski-like operators at second order in derivatives: where G µν is the Einstein tensor,R µανβ the dual Riemann tensor, while Λ is the strong coupling scale. These operators lead to consistent, second order equations of motion. They control the Lorentz-violating, timedependent scalar background profile that causes cosmic acceleration, as found in eq. (3) in the scalar example. They represent the first terms in the series of higher derivative operators, which will determine the Lorentz invariant, UV completion of the theory 2 . However, at lower energies below the scale Λ, the operators of eq. (7), together with the scalar background profile, lead to a modification of the speed of GWs.
where β = φ H/Λ 3 . Note that causality considerations do not necessarily demand δc > 0, as this statement is frame dependent [52,53]. For simplicity we have considered an example where the EM is minimally coupled to gravity and does not mix with the DE field; but in principle we could imagine generic situations where both EM and GWs differ from perfect luminality at low energies, and recover luminality at higher frequencies within the LVK band. Note however that in practice, EM waves are observed at much higher frequencies than GW ones.
In what follows we shall consider the cutoff Λ of the EFT to lie precisely between LISA and LVK scales, so that at LISA scales, the speed of GWs can be taken to be c T (k Λ) ∼ c EM − δc, while at LVK scales high-energy effects have to accounted for. At higher energies, higher frequency corrections are relevant and the expansion (8) is no longer valid. Instead a microscopic description takes over, for which the details of the spontaneous breaking of Lorentz invariance at low energy are expected to be irrelevant, and the speed of GWs is expected to be luminal (just like high-frequency light is unaffected by the medium in which it propagates). The transition between the low-frequency regime (for instance between (8)) and the high-energy one are parameterized by the details of the UV completion and can in principle be arbitrarily sharp.
Any modification of the dispersion relation can be already strongly constrained at LVK frequency scales, as proved in many analysis [54][55][56][57][58][59]: see [49] for the most recent strongest constraints. We shall therefore take an agnostic approach and focus our attention on the extreme situation where the speed of GWs is considered to be exactly luminal at LVK scales and precisely constant (but non-luminal) at LISA scales, with a sharp transition in between.
In the specific case of cubic Galileons/Horndeski systems, where the quantity G µν of eq (6) reads G µν = G 3 (φ, (∂φ) 2 ) φg µν , or for theories beyond Horndeski, characterised for instance by G µν = F 4 (φ, (∂φ) 2 ) µαβγ να β γ ∇ α ∇ α φ∇ β ∇ β φ, a resonant decay of luminal GWs into a subluminal DE field has been pointed out [60], together with non-linear instabilities that potentially rule out such models at LVK energy scales [61]. However, these bounds depend on the precise details of the EFT and its UV completion. In the present work we focus on different choices of operators, and we additionally rely on a suitable UV completion to ensure that the speed of both GWs and the dark energy field are luminal at LVK frequency scales. While at LISA energy scales the speed of GWs still departs from luminality, at LISA frequencies the decay of gravitational waves into dark energy would be strongly suppressed, as it scales with the fifth power of the momentum. As for the nonlinear instability found in [61], it is specific to the cubic Horndeski and beyond Horndeski models; but one may worry that similar non-linear instabilities may also occur for a more generic non-minimal coupling in G µν , and thus affect our results. However, as indicated in [61], the fate of this non-linear instability ultimately depends on the UV completion of the low-energy EFT, and in this work we precisely consider situations for which the UV completion (or partial UV embedding) leads to a decoupling between the dark energy field and GWs at higher energy (as expected from a standard Lorentz-invariant high-energy completion).
Even in cases where the UV completion may be ineffective at fixing such an instability (e.g. if the latter is occurring at very low energies), it is still possible that the instability is simply due to a breakdown of the Cauchy (initial-data) problem. In fact, the instability of [61] manifests itself as the effective metric on which the scalar waves travels becoming singular, which is reminiscent of the Tricomi/Keldysh breakdown of the Cauchy problem in K-essence [44-46, 62, 63] and non-linear Proca [64] theory. The latter has been shown to be fixable by suitable choices of gauge and self-interactions [45,63], and/or by a UV completion [46,65].

III. GW DATA ANALYSIS AND WAVEFORM MODEL
Having explained the theoretical motivation behind the model we study, in this section we describe the set-up of our forecasts and the parameter sets we consider.
The data d(t) (in time domain) measured by any direct GW wave detector are expressed as a combination of signal s(t), plus detector noise n(t). The quantity s(t) can be further expanded as the product of the GW signal h ij (t) with the detector response function r ij (t), which depends on the instrument under consideration: see e.g. [66] for a review on this topic. The interferometer response to a signal depends both on GW frequency A schematic representation of the frequencydependent model of cT considered in this paper. and direction. For long-duration signals, the detector position changes over time 3 . In order to avoid such timedependence, we split the total observation time T into segments of duration ∆T , with ∆T T M , where T M is the timescale of the detector motion. In this way, the detector can be assumed to be at rest within each time segment. Under these assumptions, the data can be expressed in frequency as: where the index c runs over the number of data segments. The condition of reality imposess c (f ) =s * c (−f ), and analogously for the noise. For the sake of simplicity, in the following we drop the c index.
Assuming the noise to obey a Gaussian statistics, the log-likelihood for the noise residual can be expressed as: where θ is the vector of parameters used to characterize the signal,s th (f, θ) corresponds to the theoretical model, and (a|b) is the noise weighted inner product Forecasts on the determination of the model parameters θ can be obtained in terms of the FIM approach [67,68]. Assuming the minimum of the likelihood matches with the injected value, θ 0 , the FIM -corresponding to the second derivative of the loglikelihood evaluated in the best fit -results The signal model, and in particular the GW waveform h th (f, θ), can be quite complex depending on the desired level of accuracy. In order to avoid expensive numerical relativity computations, several phenomenological templates, with increasing level of precision, have been introduced [69][70][71][72][73][74][75][76]. Typically, these templates are based on perturbative expansions (for reviews on post-Newtonian theory and on BH perturbation theory see respectively [77] and [78]) to describe the inspiraling and ring-down phases with an intermediate merger region tuned to numerical relativity simulations. In general, the total number of parameters can be quite large, with a minimum of 10 parameters for aligned and nonprecessing spins described by a single effective parameter. In the following, and similarly to [79], we adopt a simplified model where we perform an average over the sky-localization of the source, as well as over its inclination and polarization. Moreover, we do not include the effects of the spins of the merging bodies. Under these assumptions, the parameter vector reduces to: The first five (pure GR) parameters are (the log of) the observed chirp mass of the binary (in solar mass units), (the log of) its symmetric mass ratio, the luminosity distance of the binary, the time to coalescence (in ms) and the coalescence phase. Since we are assuming a sharp transition for the GW speed occurring outside both of LISA and LVK frequency bands (see Fig 1), θ MG consists of a single parameter δc, controlling the height of the transition. The exact location of the transition in frequency is irrelevant, as long as it occurs between the LISA and LVK frequency bands, so we do not include it as a parameter.
The waveform model adopted in this work is the same as described in [79], and we refer the reader to this work for details. Our model systematically expresses modifications of amplitude and phase of the waveforms in terms of functions depending on the frequency-dependent GW velocity. For the analysis presented in this paper we consider the lowest order approximation for the amplitude part of the waveform and we expand the phase up to 2.5PN order. Notice though that in the present context we can also exploit the following simplifications: • The transition between the GW speeds is sharp, and we assume it to take place outside (and in between) of the sensitivity bands of either detectors. This condition simplifies considerably the expressions for the waveforms discussed in [79]. Relaxing the sharpness of the transition would only lead to more stringent constraints, since a modification of the dispersion relation would then be potentially observable in one of the frequency bands. In fact, the transition should be very sharp for ensuring that the GW speed satisfies current bounds at LVK scales: see Appendix D of [79] for a discussion of this topic.
• Since the dominant factor in constraining variations in the propagation speed comes from the comparison of measurements at the two detectors -and not from each of the two separate observations -we further simplify the waveform, by only considering the common inspiraling part of the signal. We assume here clocks at the LISA and terrestrial detectors are synchronized to high accuracy.

IV. NUMERICAL RESULTS AND THEIR INTERPRETATION
We apply the techniques outlined in the previous section to an explicit example, aimed at verifying the analytical estimate discussed in Section I. For this purpose we consider an event with characteristics similar to GW150914, i.e. with detector frame masses m 1 = 36M , m 2 = 29M and at redshift z = 0.09. The fiducial values of t c and Ψ c are chosen to be zero. For forecasting LVK constraints, we use public O3 sensitivity curves for  LIGO Livingston, LIGO Hanford, and Virgo 4 . For our LISA constraints we make use of an effective sensitivity curve [80] obtained after averaging over sky and polarisation angle.
As a first step, we provide the constraints that this event imposes on GR parameters. We assume that the minimal detectable frequency is f m 0.017 Hz, implying that we start to observe the system 4yrs before the merger takes place. A plot of the constraints obtained from LVK observations, LISA, and a joint multiband observation of the event are shown in Fig. 2, while the specific constraints on parameters for each case in Fig. 2 are reported in Table I. We find that tighter constraints on the chirp mass and the mass ratio are obtained in the LISA band -likely because of the long inspiral stage taking place at low frequencies, resulting in a phase increase which strongly depends on M and η. However the LVK band provides tighter constraints on the other parameters. This is due to a higher SNR in the LVK band, SNR LVK = 55.69, compared to SNR LISA = 4.7 in the LISA band. As expected, the joint constraints with both detectors reach a tighter level than each detector individually.
We then proceed by considering models which allow for a sharp transition in c T to occur at a frequency between the LVK and LISA sensitivity bands. For definiteness, we choose the transition frequency to be at 1 Hz. Before studying a non-GR injection, i.e. with δc T /c EM = 0, we first perform a FIM parameter estimation with a non-GR model on a pure GR injection, i.e. δc T /c EM = 0. This particular example, whose results are shown in Fig. 3, is considered in order to give a numerical validation to the analytical estimate provided in eq. (1). For this event, as we can learn from the last row of the plot showing the constraints on δc/c EM , deviations from GR are constrained to the 10 −16 level. Besides constraints on the parameter δc/c EM , Fig. 3 also shows the forecast constraints on M o , η, d L and Ψ c . As expected, given that an additional parameter is involved in the model, the constraints slightly worsen with respect to the pure GR case summarized in the last line of Table I. However, their size keeps the same order of magnitude as in the GR example.
The final case we consider in this work is a non-GR injection. In particular, we will consider the case where the height of the transition taking place at f = 1Hz is δc/c EM = 10 −15 . A plot of the forecast constraints of this event for a joint observation with the LVK detectors and LISA are shown in Fig. 4. In this plot, the values corresponding to GR (i.e. δc/c EM = 0) are shown as a dashed black line in the bottom row of plots. We find for the 1σ bound of δc/c EM the value 1.93 × 10 −16 , which means that our fiducial value δc/c EM = 10 −15 can be comfortably detected, and the GR value δc = 0 excluded at more than 5σ level. As for the pure GR injection, the forecast constraints on all the other parameters are only mildly affected by the introduction of the additional parameter.
As described in [49], a deviation of δc/c EM = 3 × 10 −15 would be sufficient to shift by 2 minutes the coalescence time t c in the LVK band. Given our FIM forecast on t c for GW150914-like events, in such a case we will be almost certain to detect the time shift with respect to the LISA prediction. This is compatible with the results in [49] obtained using a time-shift method.
Notice also that if the transition in c T (f ) occurs within the LISA band, constraints on δc can be obtained with the GW signals from binary massive black hole mergers [79]. However, if the sharp transition takes place at frequencies lower than the LISA band, the joint detection with LIGO-Virgo-KAGRA and LISA will not be able to detect deviations on c T , since δc = 0 within their sensitivity bands. Such a low-scale transition would occur in scenarios with a parametrically very low cutoff scale as could be motivated from high-energy bounds. Constraints on low-energy EFTs coming from consistency of the high energy completion (in particular improved versions of positivity bounds [81,82] and recently derived fully crossing symmetric positivity bounds [83][84][85]) indicate that the cutoff in theories related to (7) is often required to be much lower than the expected strong coupling scale Λ.
To summarise, the forecasts performed here -using a 2.5PN order waveform -indicate that fractional changes to the speed of GWs as small as 10 −15 can be cleanly detected with a single multiband GW system. Furthermore, the remaining 'standard' parameters of the system are not degraded in any significant way, allowing for reliable analysis in a multiparameter space.

V. CONCLUSIONS
The propagation of GWs provides an exquisite probe into the nature of gravity and its couplings with the dark sector (see [86][87][88][89][90][91][92][93][94][95][96][97] for other GW propagation-related work). Motivated by the late-time acceleration of the Universe, EFTs of DE and infrared EFTs of gravity have emerged as systematic ways to parameterize couplings between gravity and the fields driving the acceleration. In the regime of validity of these EFTs (or in their weakly coupled regions), the speed of GWs can typically differ from that of the EM counterpart. However at higher frequencies, the degrees of freedom driving the late-time acceleration of the Universe are expected to be insignificant, and hence their effects on the propagation of GWs are potentially undetectable by observations in the LVK band alone.
In this work we explored the possibility that the propagation of GWs undergoes a sharp transition between the LISA and the LVK frequency bands. We proved that a single SOBHB merger event detected in both bands would be sufficient to either confirm a change in speed or severely constrain any scale-dependent transition between 10 −4 and 10 4 Hz. This conclusion can be reached analytically, but we have corroborated it by a FIM analysis of future joint LISA-LVK observations. Our analysis confirms that the observation of a single event in the both the LVK and LISA bands would allow for constraining the speed of GWs within the LISA band at the 10 −15 level (at ∼ 5σ).
Such a constraint would strongly reduce the allowed region of parameter space of the EFT of DE, and further limit models of modified gravity with modified dispersion relations. Besides testing the nature of DE (and its coupling with gravity), a multiband constraint on the propagation of GWs would provide a significant new insight on the nature of gravity in the low-energy regime and its couplings (self-interactions or couplings with other light degrees of freedom). In particular it would be interesting to explore how birefringence effects highlighted in [98] would further constrain model of modified gravity.
The rate of heavy SOBHBs detectable in the LISA band remains highly uncertain at present [28]. However, the results laid out here indicate that a single multiband system will have a significant impact on the field, analogous to the first multimessenger event and settling questions raised in its wake. In this Appendix we expand on the arguments outlined in the main text, relating the difference δc = c T − c EM with the difference |t Lisa − t Earth | between the 'expected' (by LISA) and 'measured' (by a ground-based detector) times of coalescence. During a binary inspiral and merging process, the change in the orbital radius renders the GW frequency time-dependent. At leading order in a post-Newtonian (PN) expansion, we find the relation (see e.g. [79]) with f the GW frequency measured at the detector position. t ret indicates retarded time with t the time measured at detector position, D = a(t)r the distance of the merging binary (r being the comoving distance, and a(t) the scale factor) and M is the binary chirp mass. The function F depends on the GW speed c T , and is assumed to be itself a function of frequency. The structure of F depends on the scenarios and approximations considered: see for example [79] for more details and explicit examples; in the limit of GR, F = 1. In deriving equations (A1) and (A2), we work in a nearly-static regime where the explicit time-dependence of the scale factor can be neglected, as appropriate for low-redshift sources. Moreover, we assume that δc/c is extremely small, so that in first approximation c T does not influence the amplitude of the waveform.
We can re-arrange and integrate both sides of eq. (A1), and obtain 1 If c T were one, the second terms in both sides of this equation would vanish. One would find the standard GR expression for the time τ GR = t − t c to coalescence: τ GR = 5 96 π 8/3 M 5/3 f fc df f 11/3 .
Assuming standard GR formulas, LISA observations would then estimate the coalescence time t LISA . Instead, the actual detection of the coalescence process with ground-based experiments would measure the effects of a frequency-dependent c T (if present), and a coalescence time t Earth different from what predicted in GR. Taking the difference among the two cases, we find where, since coalesce occurs at ground-based detector frequencies, we assume c T (f c ) = c EM . Starting from eq. (A5), we can estimate the order of magnitude for the bounds on δc offered by multiband measurements. We expect that the uncertainties on measurements of the coalescence time are of order O(1) s. Considering an event as GW150914, occurring at a distance of order D = 410 Mpc, and assuming it is first detected in the LISA frequency band, we have a time to coalesce of order τ GR = 5 years. Hence, D/c EM ∼ 4.2 × 10 16 s, while τ GR 1.6 × 10 8 s. Making the hypothesis that the contribution within square parenthesis in eq. (A5) is at most of order one -as expected for physically motivated choices of F, see e.g.
as stated in the main text.