key: cord-0581637-thi5ubau
authors: Datta, Sayantani; Gupta, Anuradha; Kastha, Shilpa; Arun, K. G.; Sathyaprakash, B. S.
title: Tests of general relativity using multiband observations of intermediate mass binary black hole mergers
date: 2020-06-22
journal: nan
DOI: nan
sha: 4b880907affb13d76c98992e16cb6748ca914076
doc_id: 581637
cord_uid: thi5ubau

Observation of gravitational waves (GWs) in two different frequency bands is referred to as {it multiband GW astronomy}. With the planned Laser Interferometric Space Antenna (LISA) operating in the $10^{-4}-0.1$ Hz range, and third generation (3G) ground-based detectors such as the Cosmic Explorer (CE) and Einstein Telescope (ET), operating in the $1$--$10^4$ Hz range, multiband GW astronomy could be a reality in about a decade. In this paper we present the potential of multiband observations of intermediate mass binary black holes (IMBBHs) of component masses ${sim}10^2$--$10^3M_{odot}$ to test general relativity (GR). We show that mutiband observations of IMBBHs would permit multiparameter tests of GR---tests where more than one post-Newtonian (PN) coefficient is simultaneously measured yielding more rigorous constraints on possible modifications to GR. We also find that the improvement due to multibanding can often be much larger than the best of the bounds from either of the two observatories. The origin of this result, as we shall demonstrate, can be traced to the lifting of degeneracies among the various parameters when the information from LISA and 3G are taken together. We obtain the best multiband bounds for an IMBBH with a total redshifted mass of $200M_{odot}$ and a mass ratio of 2. For single-parameter tests, this system at 1 Gpc would allow us to constrain the deviations on all the PN coefficients to below 10% and derive simultaneous bounds on the first seven PN coefficients to below 50% (with low spins).

Einstein's general relativity (GR) has been subjected to a plethora of tests performed both in the laboratory as well as using astrophysical observations [1] . The theory has so far been consistent with each of these tests (see Refs. [2] [3] [4] [5] [6] for an overview of various astrophysical tests of GR). The first observation of gravitational waves (GWs) from the binary black hole (BBH) merger GW150914 [7] and several others [8] [9] [10] [11] [12] [13] [14] during the first and second observing runs, have permitted tests of GR in a regime of strong gravity and high curvature which were elusive till date [15] . The binary neutron star merger GW170817 [16] further facilitated tests of strong-field gravity for non-vacuum spacetimes [17] . The observation of electromagnetic emission associated with this event helped in deriving stringent constraints on the speed of the GWs [18] . All these tests have placed tighter constraints on possible deviations from GR [19] , while ruling out modified theories of gravity invoked to explain dark energy [20] [21] [22] [23] [24] .

Parametrised tests of GR [25] [26] [27] [28] [29] [30] , are among the pioneering tests of the theory performed with GW data. These tests make the best use of the structure of the GW phase evolution from the post-Newtonian (PN) approximation to GR [31] . In the PN approximation, the phase evolution of the GW signal can be expanded as a power series in v and log v, where v denotes the velocity parameter describing the orbital motion of the binary. The different PN orders (corresponding to different powers of v) capture the diverse physics and various nonlinear effects * sdatta94@cmi.ac.in underlying the compact binary dynamics. Hence looking for deviations in the PN coefficients is equivalent to constraining different physics that goes into them [32, 33] . In this framework, deviations from GR are parametrized via deformation in the phasing formula at different PN orders [29, 30] whose values are put to test using the GW data. As these deformation parameters take the value zero in GR, this null test is devised to derive constraints on them at a fixed credible level.

The parametrized tests of GR branch out into several subclasses depending on the number of PN deformation parameters that are simultaneously estimated from the data. Ideally, one aims to constrain all or several of the PN deformation parameters simultaneously using the GW data [25] . This will be referred to as multiparameter tests in this paper. One may wish to further classify these multiparameter tests into two classes, depending on whether the block of PN parameters that are tested start from the lowest PN order (in the ascending order) or from the highest PN order (in the descending order). The former would make sense in terms of verifying the predictions of GR at different PN orders with increasing levels of complexities in the nonlinear interactions. The latter perspective, starting from the highest PN order and goes in the decreasing order, would be expected from modified theories such as an effective field theory where modifications to GR would start at a particular PN order and all orders above that [34, 35] . There could be other possible combinations of the PN deformation parameters that may be tested simultaneously, but we consider only these two classes of the multiparameter tests in this paper as they are the most general ones. These classes of tests, though more rigorous, yield weaker bounds, compared to single-parameter tests, due to the large correlations of the deformation parameters among themselves as well as with the intrinsic parameters of the binary such as its masses and spins [25, 36] .

Hence, one considers a somewhat less rigorous set of tests where only one of the many PN deformation parameters is chosen at a time as a test parameter [26, 29, 30] . This is less rigorous in the sense that, a modification to the phasing formula from a non-GR theory is likely to occur at more than one PN order. This aspect is not accounted for in the formulation of the single-parameter tests. This drawback is partially compensated by performing a set of tests by varying the PN deformation parameter systematically from 0PN till 3.5PN one-by-one (see [27] for a detailed discussion). Hence, one or more of these tests would potentially detect a deviation if the underlying theory of gravity is not GR, though a deviation seen at a particular PN order in this test does not necessarily mean the breakdown of GR occurs only at that particular order. Given the sensitivity of current generation of GW detectors, this is the method presently being employed in the analysis of the LIGO and Virgo data and will be referred to as singleparameter tests in this paper. The current constraints from the single-parameter tests, at 90% credible level, from the BBHs observed during the first and second observing runs of LIGO and Virgo detectors are reported in [15] . However, with the next generation ground-based and space-based detectors, the sensitivity would reach to levels where the multiparameter tests would be possible [37] .

Several studies have quantified the projected bounds on these PN parameters using third generation (3G) ground-based GW experiments such as Einstein Telescope (ET) [38] and Cosmic Explorer (CE) [39] as well as space-based Laser Interferometer Space Antenna (LISA) [40] (see, e.g., [25] [26] [27] 41] ). The ground-based 3G detectors are sensitive to stellar mass BBHs (up to ∼100M ) and intermediate mass BBHs (IMBBHs) (∼10 2 − 10 3 M ) in the frequency range ∼1 − 10 4 Hz [42, 43] , while the space-based LISA mission is most sensitive to supermassive black hole mergers (∼10 5−7 M ) in the 10 −4 − 0.1 Hz band [44] . While CE/ET will have the advantage of around 10-40 fold improvement in strain sensitivity compared to the present generation detectors such as advanced LIGO and advanced Virgo, LISA will open up the possibility of using mHz band for GW astronomy. But neither of them may be able to set stringent enough bounds on the deformation parameters to rule in or rule out viable modified theories of gravity [37] . This is because the intrinsic degeneracies of the deformation parameters with masses and spins prevent making precise measurements of these parameters using either ground-based or space-based experiments alone [27] .

In the past few years, an alternative strategy to combine these two classes of observations has been proposed as a new tool to probe the strong-field dynamics [45] [46] [47] [48] [49] [50] [51] of BBHs. This is often referred to as multiband GW astronomy [52] where, using a class of sources visible in both LISA and 3G, one combines the low frequency (early dynamics of compact binaries) content in the LISA band and high-frequency content (carrying an imprint of the late time dynamics of compact binaries) in the 3G band to obtain bounds on departures from GR. Various studies on multiband parameter estimation mostly used stellar-mass BBHs like GW150914, which will have a signal-to-noise-ratio (SNR) of order unity in the LISA band but several hundreds to thousands in the 3G band. Even then, in these studies, joint observation have been argued to be able to provide bounds several orders of magnitude better than from the individual observations [47] [48] [49] [50] [51] . This huge improvement has been broadly attributed to the combination of the low-frequency sensitivity of LISA with the high-frequency response of the 3G detectors. In Ref. [51] , these generic features have been confirmed, for the first time, with a treatment of the problem within the Bayesian inference framework.

In this work, we take this paradigm forward by carrying out an extensive study of the effect of multiband observations of IMBBHs, as opposed to stellar mass BBHs, using LISA and CE/ET, with total source frame masses of the binaries varying between 100 − 550M 1 . Astrophysically IMBHs can have masses as high as 10 4 M [53] and would be excellent sources of GWs for these detectors, although they are yet to be unravelled by observations. Therefore, the detection and understanding the formation of IMBBHs are among the top science priorities of present and future astronomical telescopes [54, 55] . The multiband detectability of IMBBHs is discussed in detail in Ref. [56] . Further, possible implications for multibanding of IMBBHs for parameter estimation and tests of GR are highlighted in [57] . In this paper we present a detailed study of the implications of multiband observations of IMBBHs by LISA and CE/ET detectors in terms of tests of GR.

Unlike the stellar-mass BBHs, the IMBBHs will have SNRs of the order of tens in the LISA band (as opposed to order of unity for stellar-mass BBHs) while in the CE/ET band they still have SNRs of the order of hundreds. This significantly helps the process of multibanding thereby providing us precise measurements of the PN deformation parameters. Further, our detailed study reveals that the dramatic improvements from multibanding is due to large scale cancellations of correlations among different parameters in the problem, due to the mutually complementarity of the two experiments or frequency bands. We find that an intermediate mass BBH with a redshifted mass of 200M and a mass ratio of 2, at a luminosity distance of 1 Gpc will yield the best multiband, multiparameter bounds on the PN deformation parameters. Considering such a system, we find that single-parameter tests can constrain first three PN deformation parameters to an accuracy ∼ 0.1% and the rest to below 10%. Multiparameter tests up to 7-parameter case can be performed with the above system, where the first two PN deformation parameters are bounded to below 0.5% and the rest to below 50% (with low spins).

The remainder of the paper is organized as follows: In Sec. II we discuss the basic concepts in combining the information from the two frequency bands (LISA and 3G), specific to the case of parametrized tests of GR. The results for the singleparameter tests with explanations of the trends seen are presented in Sec. III. Multiparameter tests are discussed in Sec. IV. Lastly, some of the caveats of the analysis are listed in Sec. V and our conclusions are provided in Sec. VI.

A. Parametrized tests of GR using IMRPhenomD waveform

The breakthrough in numerical relativity [58, 59] has enabled us to construct analytical or semi-analytical waveforms which account for the inspiral (early phase of binary evolution), merger (the late stages of the binary evolution as the two objects coalesce) and ringdown (the post-merger phase of the remnant black hole) phases [60] [61] [62] [63] [64] [65] [66] . An important subclass of them, referred to as inspiral-merger-ringdown phenomenological waveforms or IMRPhenom, are constructed starting with an ansatz about the structure of the frequency domain gravitational waveforms, which contain several free parameters that are fixed by matching with numerical relativity simulations for various mass-ratios and spins. Here we use the IMRPhenomD waveform model [62] of the IMRPhenom family, which assumes the spins of the binary constituents to be aligned or anti-aligned with respect to the orbital angular momentum vector and hence the binary is non-precessing. The amplitude of IMRPhenomD accounts for only the leading quadrupolar(l = 2, |m| = 2) mode. Schematically a frequency domain waveform would read

where A( f ) and Φ( f ) denote the amplitude and phase of the gravitational waveform. The phase, in the inspiral regime, admits an expansion of the form [31, [67] [68] [69] [70] [71] [72] ]

where t c , φ c are kinematical parameters related to the time and phase of arrival of the signal at the detector and v ≡ (πM f ) 1/3 is the PN expansion parameter in terms of which the amplitude and phase are expressed. Furthermore, η ≡ m 1 m 2 /M 2 is the symmetric mass ratio where M ≡ m 1 + m 2 is the total mass of the binary. Note that these masses are the detector frame(or redshifted) masses after accounting for the redshift of the source and related to the source frame masses M source by M = M source (1 + z). The PN coefficients in the phasing formula, deformations in which we are interested, are denoted by φ k and φ kl .

For the present analysis, we use the amplitude A( f ) and phase Φ( f ) of the IMRPhenomD waveforms, which by construction agree with the predictions of PN theory in the inspiral part of the waveform. The inspiral part of the GW phasing is described by the PN phasing formula, in Eq. (2.2), correct to 3.5PN order (O(v 7 )). Every PN coefficient φ k and φ kl are functions of the various combinations of the intrinsic parameters of the system, such as the total mass M, symmetric mass ratio η and dimensionless spins χ 1,2 of the binary components. Any modification to these PN coefficients would potentially arise from modifications to GR, an assumption which forms the basic premise for the parametrized tests. We take the ansatz for PN expansion present in IMRPhenomD and introduce free parameters at every PN order to model non-GR modifications, rewriting the PN coefficients as,

where δφ k (k = 0, 2, 3, 4, 6, 7) and δφ kl (kl = 5l, 6l) are the fractional non-GR deformation parameters to the respective PN orders denoted by index k and kl. More specifically, 5l and 6l represent the deformations to the logarithmic terms at 2.5PN and 3PN, respectively. Since the gravitational waveform is intrinsically parametrized by M, η and χ 1,2 as discussed above, our parameter space is 15 dimensional, consisting of seven GR parameters (including the luminosity distance D L of the source) and eight non-GR deformation parameters {δφ k }:

where, we find it convenient to use the chirp mass, M c ≡ η 3/5 M instead of total mass M as one of the mass parameters. We carry out parameter estimation of the GW signal described by these parameters using the projected sensitivities of LISA and 3G detectors, details of which are discussed in the next sections.

Following the huge success of the LISA-Pathfinder [73] , which has set a benchmark for a milli Hertz GW experiment in space, the European Space Agency has selected LISA for its L3 mission [44] . The proposed mission, to be launched in 2034, has an equilateral triangular constellation of three spacecraft separated by 2.5 × 10 6 kilometers, connected by six laser links. The constalletion would orbit the Sun, trailing behind Earth's orbit and has an inclination of 20 • with respect to the ecliptic. The orbital motion of the spacecraft around Sun is important for source localization and luminosity distance estimation [74] .

For our purposes we ignore the orbital motion of LISA. This is justified as our aim is to study the error on the intrinsic parameters of the binary, including the PN deformation parameters, which are more or less uncorrelated to the extrinsic parameters (such as luminosity distance, source position and orientation) and hence have minimal effect on our estimates [75, 76] . The noise power spectral density (PSD) we employ for LISA can be found in Eq. (1) of [44] , which consists of the instrumental noise and the confusion noise due to galactic double white dwarf binaries that limits LISA sensitivity in the lower-frequency regime. The latter is accounted for only 2 years of integration time in [44] . Since we use up to five years of integration time, and as the PSD scales linearly with the observation time, we multiply the expression for S gal in [44] by 5/2. The LISA noise PSD given in [44] is averaged over the sky and polarization angles and takes into account its triangular shape [74] [75] [76] [77] as well as the two independent low frequency channels. However, to account for averaging over the inclination angle which specifies the orientation of the source with respect to the detector, it will involve an additional prefactor of √ 4/5 multiplying the amplitude of the IMRPhenomD waveform model [77] . The choice of the low frequency cut-off for IMBBHs in the LISA band will depend on the duration over which the signal would last in the LISA band. If an IMBBH with chirp mass M c is observed for a duration T obs by LISA, the low frequency cut-off may be chosen as

(2.6) The Max argument above ensures that the low-frequency cutoff is not lower than the nominal low frequency cut-off of the LISA instrument. In our analysis we take T obs to be 5 years. The upper frequency cut-off is chosen to be 0.1 Hz which is equal to the upper frequency cut-off the LISA instrument.

We consider Cosmic Explorer and Einstein Telescope as two prototypical detectors representing the sensitivities that will be achieved by the third generation ground-based detectors. Cosmic Explorer is a third generation ground-based GW detector studied in the US [39] . It is conceived to be a 40 km L-shaped detector whose science goals are reviewed in [39] . We use a noise PSD for CE given in [78] S

A similar observing facility is also envisaged in Europe called the Einstein Telescope [79] . It is a triangular shaped detector with 10 km arms and is effectively three V-shaped detectors with an opening angle of 60 • . Both CE and ET noise sensitivities are limited by gravity gradient noise in the low frequency regime. The lower cut-off frequencies for ET and CE are chosen to be 1 Hz and 5 Hz, respectively. The upper cutoff frequency is chosen such that the characteristic amplitude (2 f |h( f )|) of the GW signal is lower than that of CE/ET noise by 10% at maximum. Though theoretically the upper frequency limit is infinity, contributions from the high frequency part of the waveform that contribute negligibly are ignored. This choice of the cut-off helps in significantly improving the accuracy of the numerical analysis. Similar to the case of LISA, our aim is to study the parameter estimation problem using ET and CE from a test of GR standpoint. As we use single detector configurations, we multiply the amplitude with a prefactor 2/5 to account for the averaging over the antenna pattern functions [75, 76, 80] . Figure 1 shows the strain sensitivity (square root of the noise PSD) of LISA, CE and ET. It also shows the frequency domain characteristic amplitudes (which signifies the strength of the GW signal) of an IMBBH with a total mass of 500M at 1 Gpc and a GW150914-like stellar mass binary black hole system at 400 Mpc for comparison. As the binaries inspiral, the IMBBH spends about a few hours and the stellar mass BBH a few days before they leave the LISA band at 0.1 Hz and enter ET (CE) at 3 Hz (5 Hz). The strength of the GW signals at frequencies corresponding to their respective last stable orbits is marked in magenta. This shows that the inspiral phase of systems with source frame masses greater than 500M at a luminosity distance of 1 Gpc, will hardly be visible by CE/ET. The following section discusses this in more detail. 

Following Refs. [46, 52] , there were several works which looked into the detection [81, 82] and parameter estimation [50, 56] problems in the multiband context. In Ref. [56] the authors showed that ET and CE, which are likely to be operating during the lifetime of LISA mission, will lead to multiband detections of IMBBHs up to redshift of ∼5. A typical IMBBH with total mass ∼500M at a distance of roughly 1 Gpc will be observable for years in the LISA band when the inspiraling binary components are far apart from each other. This leads to the accumulation of considerable amount of SNR during the period of observation. In Fig. 2 , we show the SNR of IMBBHs as a function of total mass for LISA, CE, ET and for multiband observations (LISA+CE and LISA+ET). Though we use IMR-PhenomD for the SNR computation, we integrate the signal up to frequency at the last stable orbit (LSO) corresponding to the total mass M, given by

This choice of upper cut-off frequency helps in explaining several features in the later sections with regard to the parameterized test of GR. At a frequency close to 0.1 Hz the IMBBH signal will leave the LISA band and after a few hours it enters the ET band at around 1Hz and the CE band at 5 Hz as demonstrated in Fig. 1 . By this time the compact binary will be inspiraling at fairly relativistic speeds until it merges. The late-inspiral and merger-ringdown phase of the IMBBH evolution will accumulate SNR of the order 1000 in CE and ET bands leading to a firm detection.

The SNR in CE (ρ CE ) and ET (ρ ET ) bands initially increases till 200M and 450M respectively and then starts to decrease as the inspiral phase of the system lasts for shorter and shorter period of time beyond these masses. SNR in the ET band is more than that of the CE from around 350M due to a better low frequency sensitivity of ET between 1 − 5 Hz. Further, the SNR in LISA band (ρ LISA ) steadily increases and matches with the SNR in CE band roughly at 550M . The multiband SNR is defined as ρ MB = ρ 2 GB + ρ 2 LISA where, GB denotes a ground-based detector, either CE or ET.

Before we discuss the technical aspects of multiband parameter estimation in the next section, it is important to point out how one would approach the problem of multibanding in practice. As shown in Fig. 2 , IMBBHs would have an SNR of the order of hundreds to thousands in 3G detectors like CE and ET. This would allow estimation of their mass parameters, especially chirp mass, to incredible accuracies. Using the time of arrival of the signal in the CE band, one would hence be able to do a search for the low frequency part of the corresponding signal in the LISA band. Due to the prior detection in the ground based detectors, the threshold on SNR for a given false alarm rate could be lowered leading to the detection of sources which have SNR ∼5 [37] . A confident detection of the IMBBH in the archival LISA data, following a search using the parameters measured by ground based detectors, forms the context for multiband parameter estimation. Indeed, more careful studies are required to quantify the detectability of these IMBBHs in the LISA band (see for instance [82] ), but we do not go into the details of this detection problem in the present work. 

In this section we discuss various elements of multiband parameter estimation we employ in this work.

Fisher matrix is a well known technique used to forecast the statistical uncertainties on various parameters in a parameter estimation problem, when both the signal and noise models are known. When the noise is stationary and Gaussian, in the limit of high SNR, the square root of the diagonal elements of the inverse of Fisher matrix yields a 1σ lower bound on the errors on various parameters [83] [84] [85] . In our case, we use this technique to compute the 1σ uncertainty in the estimation of PN deformation parameters that characterize modifications to GR. These errors also translate into upper limits on the values of the deformation parameters for a given detector sensitivity and increase in the sensitivity would lead to tighter upper limits. As discussed in Sec. II A, our parameter space is spanned by seven GR parameters and up to eight PN deformation parameters, depending on how many deformation parameters are simultaneously estimated in the problem. The Fisher information matrix is defined as the noise weighted inner product of the derivatives of the gravitational waveformh( f ) with respect to the parameters θ a that need to be estimated. More precisely,

where the noise weighted inner product is defined as

where S n ( f ) is the noise PSD of the detector and a( f ) and b( f ) are arbitrary functions of frequency. The lower and upper frequency cut-offs are denoted by f low and f high and depend on the frequency sensitivity bandwidth of the detectors, which we discussed in detail in Sec. II B. Fisher information matrix also allows the use of priors about the parameters provided they are in the form of Gaussian functions [86, 87] . If Γ (0) is the Gaussian prior matrix, the resultant Fisher matrix for any detector is the sum of the prior matrix and the Fisher matrix (Γ (0) mn + Γ mn ). The details of priors choices made are discussed in Sec. III.

For an event jointly detected by a space-based detector (LISA) and a ground-based detector (CE/ET), the multiband Fisher information matrix is simply the sum of the two Fisher matrices,

where Γ GB mn denotes the Fisher matrix corresponding to one of the ground based detectors, CE or ET. The variance-covariance matrix is defined by the inverse of the multiband Fisher matrix,

where the diagonal components, C mm , are the variances of θ m . The 1σ errors on the parameters θ m is, therefore, given as,

For any IMBBH event observed both in LISA and CE/ET, the errors on each of the binary parameters returned by the multiband covariance matrix are already marginalised over the rest of the parameters by the very definition of covariance matrix. However when there is a need to study the variance-covariance matrix of a subspace of the full parameter space (such as the one spanned by the PN deformation parameters that are estimated simultaneously), one can obtain the marginalized matrix by the following well-known prescription of constructing Schur complement of the Fisher matrix. The Schur complement of a p dimensional Fisher matrixΓ p×p is given by [88] 

where Γ p×p is the Fisher matrix block corresponding to θ p parameters that are of our interest and for which we want to study the variance-covariance matrix. Γ q×q is the Fisher matrix for θ q parameters that we want to marginalise over. Γ p×q is a matrix with cross terms between θ p and θ q parameters. Before we conclude this section, we wish to clarify a subtle issue. As the prior matrix is added to both LISA and CE Fisher matrices, one may suspect over-counting as the prior matrix is featured twice in the multiband Fisher matrix. But we would like to stress that from the parameter estimation view point, this simply reflects our assumption that parameter estimation with LISA and CE/ET are two independent experiments whose outcomes are combined to gain greater insights about the dynamics of IMBBHs.

In this section, we present the results of our analysis in detail. The first set of results are for the single-parameter tests of GR where only one of the eight PN deformation parameters is estimated along with the GR parameters. These are the first estimates of the projected multiband bounds from singleparameter tests of GR with IMBBHs and complements the earlier works for stellar mass BBHs [48, 49] . For the first time, we also provide a concrete explanation, beyond the intuitive arguments, on why the multibanding improves the tests of GR.

Our results from single-parameter tests are presented in Fig. 3 where besides the GR parameters one of the eight PN deformation parameters are estimated at a time. The bounds on the various PN deformation parameters as a function of the total mass of the IMBBH using only LISA, only CE, only ET and multiband (LISA+CE and LISA+ET) observations are shown in the figure. For convenience and to facilitate comparison between systems of different masses, the luminosity distance is fixed at 1 Gpc, the mass ratio is 2 and the component spins are χ 1 = 0.2, χ 2 = 0.1. In our analysis, we employ Gaussian priors on the component spins (mean 0, variance 0.5) and φ c (mean 0, variance π). The prior on spin is a way to remain conservative about the spin magnitudes of IMBBHs while that on φ c is meant to improve the conditioning of the Fisher matrices. The maximum source frame mass for IMBBHs is considered to be 550M because beyond this mass the signal in the CE band has negligible inspiral and hence the parametrized test would not be meaningful.

First let us examine the qualitative features. In Fig. 3 , the top panel shows bounds on deformations at the lower PN orders (0PN, 1PN, 1.5PN and 2PN ) and the lower panel depicts the bounds on the higher PN coefficients (2.5PN log, 3PN, 3PN log and 3.5PN). As expected, the lower PN order coefficients are much better constrained than the higher order ones due to their dominant contribution to the dynamics of the binary.

The bounds obtained from CE initially decrease with increasing total mass but starts to increase slowly after ∼200M following an inverse of SNR trend expected from the Fisher matrix. The bounds obtained with ET stay comparable throughout the mass range with the existence of a weak minimum around the loudest system in ET band which is ∼ 350M as shown in Fig. 2 . The errors from LISA almost monotonically decrease with increasing mass following the inverse of the SNR trend. Furthermore, one also observes cross-overs between the LISA, CE, and/or ET curves for most of the deformation parameters, which is simply an imprint of the similar cross-overs seen in the SNR curves in Fig. 2 . The total mass at which this cross-over happens is different for the SNR and the deformation parameters as the latter have more complicated noise moments [87] (powers of frequency weighted by the noise PSD) that constitute the Fisher matrix. The multiband bounds with LISA+ET are comparable to that of LISA+CE specially for the lower PN deformation parameters. Some improvements are observed in the multiband bounds on the two highest PN orders (3PN log and 3.5PN) when ET is used instead of CE.

We notice that bounds onδ φ 0 ,δ φ 2 andδ φ 3 improve the most due to multiband observations for the entire IMBBH total mass range. LISA alone helps in constraining lower PN orders more than the higher PN orders as it collects more information in the low frequency regime. However CE and ET have more information than LISA on the high frequency PN orders and are able to constrain all of them to ≤ O(1) accuracy. When we combine the Fisher matrices for LISA and CE (ET), the bounds on δφ 0 , δφ 2 and δφ 3 improve approximately by a factor of 20 to 70 (10 to 30) for masses > 300M . However, the multiband bounds on higher PN deformation parameters, δφ 4 , δφ 5l , δφ 6 , δφ 6l , δφ 7 mostly follow CE/ET except at higher masses where multiband observations improve the bounds by a factor of 5 to 10 due to the high SNR in the LISA band. These substantial improvements in bounds, specifically at the higher mass regimes, may come as a bit of surprise as the multiband bounds are factor of tens better than the best bounds obtained from either of the detectors. We devote the next subsection to explain this interesting result.

We also find that we get best multiband bounds on PN deformation parameters for equal mass case, q = 1 and the bounds worsen with increasing mass asymmetry, though not drastically. The multiband bounds also improve with increasing dimen- sionless spin magnitudes, particularly for higher PN orders.

We conclude this subsection with some quantitative statements that can be read off from Fig. 3 . The best bounds due to multibanding are obtained on 0PN and 1PN phase deformation parameters that are roughly measured to ≤ O(10 −3 ) accuracy. The rest of the parameters can be estimated to roughly between O(10 −3 ) to O(10 −1 ) accuracy for IMBBHs of total mass ranging from 100 − 550M at 1 Gpc.

As we next discuss the effect of multibanding, we consider only CE as a representative of 3G ground-based detectors. This is because the differences between ET and CE are small enough and it would not make any difference to our conclusions. Figure 3 shows that the multiband observations can provide huge improvement in the bounds of δφ 0 , δφ 2 , and δφ 3 inspite of the low SNRs in LISA band. For instance, bounds on δφ 0 from LISA and CE at 500M are O(10 −2 ) as compared to the joint multiband bound which is O(10 −4 ). This two orders of magnitude improvement may seem surprising at first. However, our investigations reveal that this feature is due to the cancellation of several off diagonal terms (which correspond to degeneracies in the parameter space) when we add the Fisher matrices of LISA and CE to obtain the multiband Fisher. Due to this cancellation the inverse of this combined Fisher matrix results in errors that are significantly smaller than the ones from LISA or CE alone.

Owing to the difficulties in representing higher dimensional matrices pictorially, we focus on selected two-dimensional subspaces, which are highly correlated, and the corresponding ellipses to understand the effect of multibanding. Consider an IMBBH system of total mass 500 M with spins χ 1 = 0.2, χ 2 = 0.1 at a distance of 1 Gpc. It is intuitive to relate the area of the two dimensional ellipses for a particular detector to its ability to simultaneously measure the two parameters. The smaller the area, the better is the measurement. Likewise, the orientation or the tilt of the ellipses tells us about the sign of the correlation between the two parameters: positive tilted ellipses (whose semi-major axis subtends an angle less than 90 degrees) refer to positive correlation between the two parameters and negatively titled ellipses indicate negative correlation between the two parameters.

Let us consider δφ 2 and δφ 7 estimates for the demonstration, as multiband improvement is the highest for the former and the lowest for the latter. Figure 4 shows the 1σ confidence ellipses of δφ 2 (top panel) and δφ 7 (bottom panel) with ln M c (left panel) and η (right panel). These two-dimensional ellipses are obtained by marginalizing over the remaining parameters following the prescription in Eq. (2.13).

The marginalized two-dimensional 1σ confidence ellipses in the ln M c − δφ 2 and η − δφ 2 plane for LISA, CE and multiband (LISA+CE) are shown in the top panel of Fig. 4 . The area of the contour from CE is much larger than and a tilt that is orthogonal to the one from LISA, implying that LISA can measure both parameters better than CE and the correlation between these two parameters is positive for LISA and negative for CE. The multiband two-dimensional confidence ellipse is in black, whose zoomed-in version is shown in the inset. The area of this ellipse is substantially smaller than that of the LISA and CE. The origin of this feature may be traced to the reduction or cancellations of several off-diagonal terms when the multiband Fisher matrix is constructed. Since the marginalization involves cross-terms in the covariances matrix, the two-dimensional plots are sensitive to such cancellations. Table I provides the areas of the two-dimensional ellipses for various parameter combinations. This is a quantitative demonstration of the complementarity of the two frequency bands very often invoked in the literature to explain the improvements due to multibanding [57] . The top-right panel of Fig. 4 shows the confidence ellipse for η − δφ 2 for LISA and CE. Unlike the ln M c − δφ 2 case, though LISA seems to be able to measure the parameters better, CE's ability to measure these parameters are comparable though slightly worse (see Table I for the estimates of the area). When combined with the opposite tilts of the two ellipses, the multiband estimates help improve the parameter estimation and hence the bounds.

Similar two dimensional ellipses are presented for ln M c − δφ 7 and η − δφ 7 in the bottom panel of Fig. 4 . Here, as one may read off from Table I , the ln M c − δφ 7 ellipses have the smallest area for LISA, though the two areas are not very different as it was for ln M c − δφ 2 . This simply tells us that for higher PN deformation parameters, CE with its high frequency sensitivity, has much better ability to lift the degeneracy. The bottom-right panel clearly depicts that CE is much more efficient in breaking the degeneracy between δφ 7 and η than LISA which again is due to the high frequency sensitivity of CE close to the merger of the IMBBH, where higher PN order dynamics becomes important.

To conclude this discussion, we have explicitly shown that the tremendous improvements due to the multiband observation is a direct consequence of the cancellation of various degeneracies when we combine the information from LISA and CE.

So far we have presented the bounds on the PN deformation parameters when only one of them is estimated at a time. There are eight PN coefficients up to 3.5PN order which give us eight single-parameter tests of GR. As argued earlier, the most general test of GR we can carry out is the one where all of the PN deformation parameters are simultaneously measured. Due to the inherent degeneracies between the deformation parameters and the GR parameters, this test is not feasible using GW observations in a single frequency band as argued in [25] .

We now study how a set of PN deformation parameters may be measured simultaneously by combining observation of the same signal in LISA and CE bands. We consider two different types of multiparameter tests of GR. In the first type, we increase the number of parameters that occur at different PN orders starting from the lowest (0PN) order. The second type of tests measure the deviation in the PN coefficients starting from the highest (3.5PN) order and going to the lower PN orders. Hence this set of tests would quantify our ability to constrain possible deviations for those theories which predict the last n coefficients to differ from GR. For instance, n = 3 would be a test where we simultaneously measure the last three PN parameters in the phasing formula. Below we provide the projected bounds on a prototypical IMBBH system with a total mass of 200M for these two subclasses of multiparameter tests.

From Fig. 3 it is evident that in the case of single-parameter tests, the total mass at which the joint bounds on most of the parameters are minimum is around 200M . This can be understood from two features which have already been discussed: (i) the trends in the SNR as a function of total mass, and (ii) higher PN coefficients play a dominant role in the late time dynamics, close to the merger, which falls in the CE band. The effectiveness of a multiparameter test of GR, or for that matter any test of GR, depends on the optimization of these two effects, which leads to a sweet spot for these tests. In our case, this happens to be around a total redshifted mass of 200M . The errors on different deformation parameters in Fig. 3 show slightly different minima, hence there are other values around 200 M that are equally good for these tests. As this choice would have a negligible impact on the conclusions we draw, we stick to a total mass of 200 M with two different choices of component spins to show the projected bounds using multiparameter tests.

A. Bounds from the lower order PN side Figure 5 shows the bounds on deformation parameters obtained from the various n-parameter tests starting from 0PN, for an IMBBH with a total redshifted mass of 200M and a mass ratio of 2 with two different spin configurations: χ 1 = 0.2, χ 1 = 0.1 (top panel) and χ 1 = 0.8, χ 2 = 0.7 (bottom panel). The binary is assumed to be at 1 Gpc. As there are seven GR parameters, for each test, we invert a Fisher matrix of dimension 7 + n to obtain the corresponding errors. We find that simultaneous measurement of only seven of the eight PN deformation parameters is possible for this binary configuration if we require the errors on all the PN deformation parameters to be less than or equal to unity. Comparing the top and the bottom panels, it is evident that the increase in spin magnitudes have varied effect on the estimation of PN coefficients, sometimes improving and other times worsening the error bounds. The lower order PN coefficients are also largely unaffected by spin magnitudes. This is a reflection of the fact that spin effects are higher order effects (starting at 1.5PN order) and hence spin dynamics plays a dominant role in the late time dynamics where, again, CE sensitivity has an important role, leading to improved bounds of the higher order PN deformation parameters.

For high values of spin magnitudes, single-parameter tests on the Newtonian and 1PN coefficients would yield constrains of O(10 −4 ) while all other parameters, except the 2PN and 3PN logarthmic ones, will be bounded to O(10 −3 ). The worst bounds are for the 2PN and 3PN logarithmic terms which are of O(10 −2 ). This precision is unprecedented compared to what LISA and CE would be able to do for supermassive and stellar mass BBHs, respectively, which they are best sensitive to [37] . By increasing the number of parameters that are simultaneously measured, the bounds on all the lower order PN deformation parameters worsen due to the degeneracies present in the waveforms. The addition of the δφ 5l , δφ 6 , and δφ 6l deformation parameters have negligible effect on the bounds on the lower order PN deformation parameters as the correlations of these lower order parameters with the higher order ones are rather weak. However, adding the 3.5PN parameter significantly worsens the bounds on all deformation parameters above 2PN, making the errors go above unity and hence are not shown. This trend is a consequence of the superior ability of LISA to measure the lower order PN coefficients and CE to measure the higher order ones. As we keep adding higher order PN coefficients, the bounds on the lower order ones, which benefit mostly from LISA, are unaffected. However, when we add more and more higher order PN parameters, such as 3.5PN, CE's ability to simultaneously measure them diminishes, leading to an overall worsening of higher order PN deformation parameters' measurement. Despite this, it is impressive to note that, with multibanding, a seven parameter test of GR can yield bounds ≤ O(1) for all seven PN deformation parameters.

B. Bounds from the higher order PN side Figure 6 shows the second set of results related to the bounds on the highest n PN deformation parameters. For instance, one can constrain any modified theory of gravity which predicts de- viations from GR starting at 2PN order (five parameter bounds denoted by pentagons) with a precision O(10 −1 ). Though CE is mostly sensitive to the higher order PN coefficients, it does not have enough ability to break the degeneracy between two consecutive PN coefficients which are strongly degenerate. In this case, LISA is also able to offer little help as its role is limited to estimate M c and η very well and break their degeneracies with higher PN coefficients. As these correlations are rather weak, LISA does not help much for this class of tests.

In this section we discuss some of the caveats of our analysis. Uncertainties about the IMBBH population: An important caveat of the results presented here is the uncertainty in the merger rates of IMBBHs. As LIGO/Virgo detectors are yet to see an IMBBH, we have to rely on the state of the art upper limits on their merger rates reported from the first two observing runs [89] . This seems to suggest that IMBBHs with a total mass 200M may have a merger rate lower than ∼ 0.2 Gpc −3 yr −1 . Since the multiband distance reach for this system is ∼ 1Gpc, this would mean an upper limit on the detection rate of, roughly, one event per year. As we go to higher masses, the upper limits are less stringent. However, as shown earlier, the multiparameter test may become increasingly difficult to carry out for a total mass ≥ 550M . Hence a clear picture about masses of IMBBHs for which this test would perform well will be clearer only in the future with more stringent upper limits or perhaps a detection. The uncertainties in the merger rate and lack of a mass distribution model for IMBBHs also prevents us from assessing how well can a population of IMBBHs improve the test. Use of Fisher matrix based parameter estimation: We have relied on the ability of Fisher information matrix approach to predict the precision with which the PN deformation parameters can be estimated using LISA and CE/ET. The Fisher matrix based approach is expected to be reliable in the limit of high SNR [86, 90, 91] . However, the projected bounds on non-GR parameters from the Fisher matrices for GW150914 and GW151226 have shown good agreements with the results from Bayesian inference [19] reinforcing the utility of Fisher matrix to obtain order of magnitude estimates of the errors. In our case, the SNRs in the CE band are of the order of hundreds to thousands and hence is well within the domain of applicability of Fisher matrix. However, the LISA SNRs are of order ≥ 10 which theoretically falls only marginally within the domain of applicability of this method. Hence our LISA-only results are prone to have uncertainties which need to be quantified using numerical sampling techniques such as Markov chain Monte Carlo [92] or Nested Sampling [93] . A recent work [51] has paved the way for more work in this direction.

Neglect of precession and subdominant modes in the gravitational waveforms: The bounds reported in this paper were obtained using the IMRPhenomD waveform model, which models a non-precessing black hole binary. This model does not account for effects such as spin-induced orbital precession [94] and subdominant modes of GW signal [71, 72, 95, 96] . The incorporation of precession [63] and subdominant modes [64] bring in characteristic modulations to the phase and amplitude of the waveform and hence is believed to be more informative in improving the overall parameter estimation (see, for instance, [27, [97] [98] [99] [100] ). Therefore, one would expect our bounds to improve with the incorporation of these effects, however, this is outside the scope of the paper.

Neglect of LISA's orbital motion: Our model for LISA noise PSD does not account for its orbital motion. As these orbital modulations have negligible impact on the estimation of intrinsic parameters of the binary [75, 76] , our estimates are unlikely to be affected much by this assumption.

We have discussed in detail the possibility of multiband observation of IMBBH systems using 3G ground-based detectors and the space-based LISA. It is shown that observation of IMBBHs would be an excellent new class of sources for tests of the strong-field dynamics. Besides the single-parameter tests of GR, IMBBHs would facilitate multiparameter tests, which simultaneously measure more than one PN deformation parameter. The addition of information from LISA and CE or ET leads to significant improvements due to massive cancellations of the off-diagonal entries of the Fisher matrix signifying how multibanding helps break the degeneracies between various parameters in the gravitational waveform. We have discussed how the projected bounds would vary as a function of total mass of the system and find that IMBBHs with a total redshifted mass between 200 − 400M would be the sweet spot for the test.

EMR/2016/005594 and a grant from the Infosys Foundation. This document has LIGO preprint number P2000209. We thank the frontline workers combating the CoVID-19 pandemic without whose support this work would not have been possible.

LIGO Scientific)

relativistic Objects in Compact Binaries: From Birth to Coalescence Editor: Colpi et

Mathematical methods in statistics

International Series of Monographs in Electronics and Instrumentation

The Schur Complement and Symmetric Positive Semidefinite (and Definite) Matrices

We thank Emanuele Berti, Ssohrab Borhanian, Arnab Dhani, Muhammed Saleem and, especially, Ajit Mehta for useful discussions and comments. B.S.S. is supported in part by NSF grants PHY-1836779, AST-1716394 and AST-1708146. S.D. and K.G.A. also acknowledges partial support by the Swarnajayanti grant DST/SJF/PSA-01/2017-18 of DST-India. K.G.A. also acknowledges partial support by and SERB grant