{'fraction_non_alphanumeric': 0.03921448672458094, 'fraction_numerical': 0.017098984461030222, 'mean_word_length': 4.726572079173235, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bounding Alternative Theories of Gravity with Multi-Band GW Observations

We study the constraints on alternative theories of gravity that can be determined by multi-band observations of gravitational wave signals emitted from binary black hole coalescences. We focus on three types of General Relativity modifications induced by a generalised Brans-Dicke theory, and two classes of quadratic gravity, Einstein-dilaton-Gauss-Bonnet and dynamical Chern-Simons. Considering a network of space and ground-based detectors, supplied by a population of spinning binaries black hole, we show how the multi-band analysis improves the existing bounds on the theory's parameters by several orders of magnitude, for both pre- and post-Newtonian deviations. Our results also show the fundamental role played by an interferometer in the frequency range between LISA and advanced detectors, in constraining possible deviations from General Relativity.


I. INTRODUCTION
Gravitational waves (GW) emitted by the coalescence of black hole (BH) binaries are among the cleanest and most valuable tools to investigate the features of gravity and to test the predictions of General Relativity (GR) in the highly-relativistic, strong-field regime [1][2][3]. GW signals observed so far have been deeply analyzed by the LIGO/Virgo Collaboration [4][5][6] to find possible deviations from GR. Recently, a series of tests has been performed on the joint datasets collected by advanced detectors, showing no sign of inconsistency with Einstein's theory [7]. Future binary black hole (BBH) observations are expected to lower statistical errors on the source's parameters, while numerical and semianalytical studies promise to reduce the systematics in the waveform's modeling [2,8,9]. Moreover, a new family of GW interferometers will complement the existing detectors with both ground-based and space facilities.
By the mid-2030s, the satellite LISA [10][11][12][13] will be ready to start an observational campaign aimed to detect GW signals emitted between 10 −4 and 10 −1 Hz. With a sensitivity curve being designed to follow the evolution of supermassive black holes, either in symmetric or extreme mass ratio binaries, LISA will also provide a complementary window, in the low frequency band, for stellar-mass black holes which are among the primary targets of ground-based interferometers. The relevance of space and terrestrial joint detections has been recently investigated [14], showing the improvement in the measurements of the source's parameters, which is crucial, as an example, for a precise binary localisation. Had LISA been operating during the first LIGO/Virgo observations, this would have allowed us to know with great accuracy the moment at which the signal would have entered the bandwidth on Earth, boosting the quality of data analysis [14,15].
Multiband observations also represent a powerful approach to probe theoretical foundations of General Rel-ativity and test gravity modifications [1] in different dynamical regimes. Stellar-mass BBH observed by both space-and ground-based interferometers evolve through a wide range of frequencies emitting GW in the millihertz LISA's band for years (the early inspiral phase), before chirping at high frequencies and producing a short signal (the late inspiral and merger phases) in the LIGO/Virgo band, around approximately 100 Hz. A joint detection would help to constrain the source's parameters, which may be dominant either at low or at high frequencies, and therefore would be measured with different accuracy by space and terrestrial detectors. For example, a pre-Newtonian effect, i.e., a non-GR correction which modifies the waveform before the leading quadrupolar order, plays a major role in the early evolution of a binary system. Strong bounds on this modifications can be placed by observing double pulsars in the electromagnetic channel [16][17][18]. For binary black holes, this effect would be more dominant in the low-frequency spectrum of LISA, leading to very different constraints with respect to a possible LIGO/Virgo detection [19].
Besides LISA and current advanced interferometers, the GW discovery is pushing the development of new facilities. KAGRA is close to completion [20,21], while huge improvements in the sensitivity will be given by third generation detectors, as the Einstein Telescope [22] or the Cosmic Explorer [23,24]. Finally, new space satellites are under active study, as the Japanese B-DECIGO [25], which is designed to bridge the gap between LISA and ground-based interferometers, and promise to observe both BH and neutron star binaries with exquisite precision 1 [27].
This network of interferometers actually represents a wideband detector able to measure with pinpoint ac-curacy the parameters of binary sources, and to test gravity throughout the full orbital evolution up to the merger phase. In this paper, we use such global web of interferometers to explore how multiband GW detections can constrain the fundamental parameters of three different alternative theories of gravity: a generalized model of the Fierz-Jordan-Brans-Dicke theory [19,28] and two types of quadratic gravity, namely Einsteindilaton-Gauss-Bonnet [29,30] and dynamical Chern-Simons [31][32][33]. In such theories, the gravitational interaction is mediated by an extra scalar field coupled with terms proportional to the curvature. The presence of the scalar field activates new physical mechanisms, such as the emission of dipole radiation, which increase the overall gravitational wave flux emitted by binary systems and hence change their orbital evolution. These corrections modify the GW phase at different post-Newtonian (PN) orders and then affect the signal in different frequency regimes. They represent therefore a test bed for the full potential of the multiband analysis.
We model the GR deviations using the parametrized post-Einsteinian (ppE) approach [34,35]. Similar to the low-velocity, weak-field PN expansion of the metric and matter variables [36], the ppE formalism maps modelindependent deviations from GR directly into the GW signals emitted by binary sources. The effectiveness of this framework to detect deviations produced by alternative theories using synthetic and real data has been deeply explored so far [19,[37][38][39][40][41][42][43][44][45][46][47] (we refer the reader to the review [48] and reference therein for an extensive lecture on the subject). In particular, Yunes et al. have explored the fundamental physics implications that current and future GW detections may have on a large set of modified theories of gravity [44,45]. Astrophysical bounds on the lowest order ppE coefficients have also been derived using observations of relativistic binary pulsars [42].
In this scenario, the ppE formalism represents a precious tool to perform multiband analysis which spans the GW spectrum within the frequency range [10 −4 −10 3 ] Hz, i.e., for a global network composed of space-and groundbased interferometers. Joint constraints on scalar-tensor theories of gravity, with a −1PN correction in the GW phase, have been studied by Barausse et al. assuming LISA and advanced LIGO/Virgo observations of prototype binaries [19]. This work has shown how the constraints on the non-GR parameters can benefit from multiple detectors, improving the bounds from single measurements of various orders of magnitude.
In this paper, we pursue a similar path, and we extend previous studies in order to (i) broaden the analysis to different classes of alternative theories of gravity that introduce distinct ppE corrections in the waveform, (ii) consider an astrophysical population of spinning stellarmass BH binaries, and (iii) take into account third generation interferometers such as the Einstein Telescope [49,50] and a second space instrument in the hertz band given by the proposed Japanese B-DECIGO [27]. We derive the distribution of statistical errors on theory's parameters, computing the projected bounds for different detector's configurations. We show how space and terrestrial interferometers, supplied by the population of BBHs, represent an incredible opportunity to test gravity modifications which affect GW signals on a wide range of regimes.

II. BINARY BLACK HOLE POPULATION
The astrophysical BBH population used for the injections is obtained by combining state-of-the-art population synthesis simulations with a cosmological simulation, as already described in [51,52]. In particular, the population-synthesis simulations provide information on BBH mass and delay time (i.e., the time elapsed between the formation of the progenitor stellar binary and the time of the merger), while the cosmological simulation outputs provide information on the merger redshift and on the host galaxy of the BBH.
We use the publicly available Illustris-1 cosmological simulation [53][54][55] with a box of 106.5 comoving megaparsecs length and a baryonic mass resolution of 1.26 × 10 6 M . The population-synthesis simulations were run with mobse, which includes state-of-the-art prescriptions for stellar winds and black hole formation [56][57][58]. In our model, the BBH mass and merger rate strongly depend on the metallicity of progenitor stars; metal-poor (Z ≤ 0.002) massive stars are more efficient in producing merging BBHs and form more massive merging BBHs (up to approximatively 90 M ) than metal-rich stars [57]. The mass spectrum and merger rate of BBHs predicted from mobse is fairly consistent with current constraints from GW detections [59,60]. The mobse simulation suit adopted here is named run "CC15α5" and was already presented by [57]. For more details, we refer to [52] and [57].
From the population-synthesis simulations, we generate catalogs of merging BBHs, which we plant in the cosmological simulations via a Monte Carlo algorithm (based on the star formation rate and metallicity evolution). From this procedure, we obtain a population of BBHs with astrophysically motivated mass and redshift distribution [61]. In these simulations, we consider only BBHs formed from isolated binaries, and we neglect other possible formation channels (e.g., dynamical evolution in dense stellar systems [62][63][64][65]. Spins are assigned to BBHs a posteriori, by randomly drawing the dimensionless spin magnitude χ from a Maxwellian distribution with root mean square equal to 0.1. Spin orientations are assumed to be aligned with the binary orbital plane. The spin model is just a simple toy model because the astrophysical processes which affect the spin magnitude of black holes are still largely uncertain (Bouffanais et al., in preparation). Our simplified description of spins is consistent with current GW data [60]. For this paper, we limit our sample to BBHs merging up to a luminosity distance of 1 Gpc (corresponding to redshift z ∼ 0.2), which is approximately the maximum distance at which LISA will be able to detect (rather massive) stellar BBHs [10,27]. The BBHs considered in this work span a mass ratio 1 m 1 /m 2 4. The distribution of the component masses for the BBHs analyzed is shown in Fig. 1.

III. WAVEFORM MODEL
The main goal of our analysis is to compute bounds on alternative theories of gravity using multiband detections of the GW signals emitted by the BBH population described in Sec. II. As discussed in the Introduction, we work within the ppE framework, which introduces modelindependent deviations from GR, both in the amplitude and in the phase of the gravitational waveform [34]. In this paper, we focus on the inspiral part of the signals, taking into account non-GR corrections in the waveform phase only, as deviations in the GW amplitude can be considered subdominant [19,66]. Extra information can be extracted from the other phases of a BBH coalescence, once the full inspiral-merger-ringdown emission in alternative theories is completely understood. New efforts have been currently devoted in this direction to study the binary evolution with GR modifications and fulfill the gap between analytical and numerical GW templates [67][68][69].
Within the ppE approach, the template in the frequency domain is given by: where u = πMf and M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 is the chirp mass of the binary [34,48].
where v is the orbital velocity. For the GR waveform h GR , we use the PhenomB template for non-precessing, spinning BBHs [70,71], assuming a Newtonian amplitude with average sky orientation.. The cutoff f IM is the inspiral-merger transition frequency defined in [71]. Finally, (b, β) are the ppE parameters, where b controls the nature of the non-GR deviations and β controls their magnitude. In the standard post-Newtonian terminology, a ppE coefficient identified by a given b corresponds to a (5 + 3b)/2 PN order term. With a specific choice b, different modifications at various PN orders, including negative pre-Newtonian effects [2], can be studied separately. In general, corrections with negative (positive) PN orders have a major impact on the low-(high-)frequency part of the GW signal. Such corrections are not completely unconstrained, since bounds in both the gravitational and the electromagnetic spectrum do exist [4,42,44].

A. Non-GR modifications
We provide here a brief description of the three alternative theories of gravity considered in the paper. The non-GR modifications that such theories bring to the GW phase cover a wide range of post-Newtonian orders. Therefore, they represent a good set of candidates to explore the feasibility of the multiband analysis, and to understand the complementarity of space-borne and ground-based interferometers to constrain a specific correction. We also focus on effects which generally belong to modifications of the GW generation mechanism [48], which are active in the source's near zone, leading to changes in the binary equations of motion. The three alternatives considered fit the general class of scalar-tensor theories of gravity. These are among the most natural modifications of GR, in which the gravity sector is nonminimally coupled with an extra scalar field. The latter affects both the binary orbital evolution and the gravitational wave flux and induces the emission of a dipole radiation, which is forbidden in GR [72]. In these theories, the Klein-Gordon equation for the scalar field is in general sourced by two terms, φ ∝ S matter 1 + S curv 2 , which are proportional to the matter's stress-energy tensor and to the curvature corrections added to the Lagrangian, respectively. The presence of both or of only one of these two components depends on the specific theory considered and determines the features of the dipole emission.
The first and simplest GR extension we focus on is a generalized version of the Fierz-Jordan-Brans-Dicke scalar-tensor theory [28], in which the scalar field is coupled to the Ricci scalar. In this theory, S curv 2 = 0, and φ is sourced by the matter content only. This implies that for a globally vacuum spacetime, for which also S matter 1 = 0, the scalar-field profile is constant. Binary black holes therefore do not lead to a dipole emission 2 . However, it is possible to have non-trivial configurations for φ for a specific choice of the boundary conditions [73]. We therefore follow the approach pursued in [19], by considering a generic parametrization which captures the effect of Brans-Dicke-like (BD-like) theories through a single coefficient B, which quantifies the magnitude of the dipole flux in the GW luminosity, dE dt = dEGR dt 1 + B M l −1 , where dE GR /dt is the quadrupolar GR component, M = m 1 + m 2 is the total mass of the binary, and l is its orbital separation. BDlike theories are described within the ppE approach by the following choice of parameters, with ν = m 1 m 2 /M 2 denoting the symmetric mass ratio, and correspond to a −1PN correction, which is therefore dominant at low frequencies with respect to the standard emission in GR.
We also consider two specific examples of quadratic gravity theories, i.e., Einstein-dilaton-Gauss-Bonnet 3 (EdGB) [29,30] and dynamical Chern-Simons (dCS) [31][32][33], in which the Einstein-Hilbert action is modified by introducing more complex corrections, quadratic in the curvature, which couple with the scalar field φ. In these theories both S matter 1 and S curv 2 are in general different from zero, and the matter-independent term leads to dipolar/quadrupolar emission in EdGB/dCS black hole binaries. Due to the nature of the non-GR terms, quadratic theories are particularly relevant in the strongfield regime of gravity. The ppE map to EdGB gravity leads to the following parameters [48]: The phase amplitude β EdGB depends on the coefficients which show a strong correlation between the BH masses and the spin parameters χ i ∈ [−1, 1] [44]. The value b EdGB = −7/3 identifies again a −1PN correction: BH binaries emit scalar dipole radiation induced by the black hole's individual monopole charges. For a given α EdGB , which is the actual parameter entering the theory's Lagrangian and has the dimensions of a squared length, the dimensionless coupling ζ EdGB = (16π/M 4 )α 2 EdGB is larger for low-mass sources. In this regard, the stellarmass population studied in this work represents the ideal arena to test EdGB theory. Finally, the last theory considered in this paper, dCS gravity, breaks two GR foundations. As for EdGB, dynamical Chern-Simons violates the strong equivalence principle. Moreover, the scalar field's equation is not invariant under parity transformation, due to the specific nature of the quadratic curvature correction, given by the Pontryagin invariant [74]. However, unlike the theories described above, dCS introduces a 2PN deformation on the GR phase, specified by the following parameters, where δ = (m 1 − m 2 )/M , χ s = (χ 1 + χ 2 )/2, χ a = (χ 1 − χ 2 )/2, and ζ dCS = (16π/M 4 )α 2 dCS [47]. In this theory the black holes feature a scalar dipole charge, which in turn induces a scalar quadrupolar emission during the orbital evolution.

IV. STATISTICAL ANALYSIS
In this paper, we consider the sensitivity of four different detectors which form a multiband GW network: i) a LIGO/Virgo second generation interferometer at design sensitivity [75]; ii) the future third generation interferometer Einstein Telescope (ET) in the so-called ET-D configuration [49,50]; iii) the Japanese space-based interferometer B-DECIGO, also planned to test the capabilities of the larger detector DECIGO [27,76]; and iv) the space-borne detector LISA, with and Optical Measurement System 4 (OMS) of 10 pm/ √ Hz at high frequencies and assuming a mission lifetime of four years [10]. The sensitivity curves of these four GW detectors are shown in Fig. 2.
The gap of approximately [0.1, 3] Hz between LISA and the ET can be filled by the decihertz band of B-DECIGO, which becomes, in this regard, crucial to construct a complete network of detectors with a seamless spectrum f ∈ [10 −4 , 10 4 ] Hz, as shown in Fig. 2.
The estimate of the source parameters, and hence of the non-GR modifications, would greatly benefit from multiband observations [19], especially for those theories which predict changes which are dominant in a specific part of the spectrum. For example, dipole radiation arising in EdGB or BD-like theories occurs at low frequencies as a pre-Newtonian correction, and it is expected to be constrained with more precision by space detectors operating at subhertz frequencies, like LISA and B-DECIGO. On the other hand, dCS gravity introduces a 2PN order effect, which becomes relevant at high frequencies, where ET plays a major role.

A. Signal-to-noise ratios and the Fisher matrix
In order to describe the mathematical tools used to compute the uncertainties on the ppE parameters, it is useful to introduce the inner product between two functions A(t) and B(t) in the waveform space, weighted on the power spectral density S n (f ) of a given detector [66], withÃ(f ) andB(f ) being their Fourier transforms. The signal-to-noise ratio (SNR) ρ of a GW signalh(f ) is then given: The bandwidth [f in , f fin ] is fixed by the specific instrument and/or GW event considered. In particular, f fin is chosen to be equal to f IM when we consider LIGO/Virgo and ET and equal to 0.1 and 100 Hz for LISA and B-DECIGO, respectively. On the other edge of the frequency spectrum, f min is limited by seismic noise for ground-based detectors, f LIGO/Virgo min = 10 Hz, f ET min = 3 Hz, while we fix f min = 0.01 Hz for B-DECIGO [27]. Finally, following [77], the initial frequency for LISA is determined by: Hz , (8) where T obs = 4 yr.
In our analysis we focus on high-SNR signals, asking that the value of ρ in each detector of the network is larger than a fixed threshold, namely ρ th = 15 for LISA [78] and ρ th = 8 for the other interferometers [79]. Figure 3 shows the SNR distributions for the population of binaries studied in this work. The requirement ρ > ρ th cuts a large fraction of events observed by LISA. This is somehow expected, as the BBH analyzed feature a stellar-mass distribution, and therefore represent the prototype target for ground-based detectors, like ET. The latter shows the largest values of ρ, followed by B-DECIGO, for which the binaries accumulate SNR due to the long sweep in the low-frequency band. Under the large-SNR hypothesis [80], we can assume that the GW parameters determined by the data analysis are Gaussian distributed around the true valuesθ, i.e., θ =θ+δθ. In this case, assuming flat priors, the posterior probabilities for the source parameters can be written as [66], where Γ ij is the Fisher information matrix, with ( · | · ) being the noise-weighted inner product (6). In this framework, the covariance matrix is simply determined by the inverse of the Fisher matrix, Σ ij = (Γ −1 ) ij , and the parameter's standard deviation given by σ i = √ Σ ii . In our analysis, we consider θ = {M, ν, χ, t c , φ c , θ ppE }, where χ is the PhenomB effective spin parameter; t c and φ c are, respectively, the time and phase at the coalescence; and θ ppE represents the ppE parameter that encodes the magnitude of the non-GR modification to the GW phase [cf. Eqs. (2), (3) and (5)].
We compute the Fisher matrix in each of the network's detectors for all the BBHs included in our catalogue, injecting the "true" values of the parametersθ in Eq. (10) under the null hypothesis that GR is the correct theory of gravity, i.e., such thatθ ppE = 0. Assuming that the GW observations are independent, it is straightforward to combine different datasets for LIGO/Virgo, ET, B-DECIGO and LISA, in order to obtain multiband constraints. In this case, where the index k runs over the N d detectors considered and Γ (k) are the associated Fisher matrices. Note that for network configurations, with LISA and one groundbased detector, eq. (11) is strictly valid only for signals which are phase connected [19]. The results presented in the next section for the multiband analysis correspond to the BBHs which satisfy this requirement.

V. RESULTS
Assuming that GR is correct, the standard deviation of each ppE parameter resulting from the Fisher matrix analysis imposes an upper bound on the corresponding modified theory of gravity at 1σ confidence level. Figure 4 shows the distributions of the 1σ errors for the parameters which characterize the three alternative theories described in Sec. III A and the BBH population of Sec. II. For each detector, we only consider sources yielding a SNR above the threshold ρ th . This requirement excludes approximately 90% of the binaries in the LISA band, while it has small impact for the other detectors. As expected, for both BD-like and EdGB, which predict a −1PN correction to the GW phase, low-frequency detectors provide the strongest constraints. For the two space detectors, we obtain: B 10 −10 and ζ EdGB 10 −11 for B-DECIGO, whereas B 10 −8 and ζ EdGB 10 −9 for LISA (1σ median values; see Table I). The most stringent constraints on B and on the EdGB coupling constant come from the Japanese detector, which compared to LISA features higher values of the SNR. In the EdGB case, it is also crucial to have very accurate measurements of the BH spin parameters [44], which are correlated with ζ EdGB (and the masses) in the GW phase through the ppE coefficient (3). In this regard, we note that B-DECIGO is able to perform measurements of χ with a relative accuracy well below 1% for all the considered BH binaries, compared against the errors obtained by LISA, which cluster around σ χ /χ ∼ 5%.
These results change when we analyze the data in terms of ground-based detectors. The left and center panels of Fig. 4 show indeed that both third and second generation interferometers would constrain BD-like and EdGB parameters several orders of magnitude less than space-borne detectors. For ET, in the best case scenario, the new bounds on ζ EdGB can be as large as 6 (7) orders of magnitude with respect to B-DECIGO (LISA). This is mainly due to the sensitivity in the frequency band of space-borne instruments with respect to the dipole correction and to the ability to constrain the BH spins with better accuracy. For BD-like corrections our results are consistent with the uncertainties on B found in [19,45] for stellar-mass binaries observed by both ground-based and space interferometers. The right plot of Fig. 4 shows the distributions of projected bounds for dCS theory. As explained in Sec. III A, in this case, the non-GR modifications lead a 2PN correction in the GW signal, which is therefore expected to contribute more at higher frequencies. The constraints coming from LISA are much looser and compatible with the results of advanced detectors, while ET is able to put more stringent bounds and to measure ζ dCS 1 (1σ median value). Note that the values obtained for LISA and LIGO/Virgo violate the small-coupling approximation, namely ζ EdGB,dCS 1, assumed to derive eq. (5). Even for dCS, the detector which provides the tightest bounds is B-DECIGO, for which we obtain a median ζ dCS 10 −1 . This result is again due to the exquisite sensitivity of B-DECIGO to the inspiral phase of stellarmass objects, which is able to constrain with good accuracy all the waveform parameters. Comparing the three panels of Fig. 4, we finally note that the constraints on GW deviations arising from BD-like and EdGB are in general several orders of magnitude stronger than those on the dCS theory. This difference is mainly related to the different PN order which characterizes the non-GR modifications in the gravitational waveform.
For both quadratic gravity theories, knowing the mass of the BBHs and their uncertainties we can propagate the errors on (ζ EdGB , ζ dCS ) into bounds on (α EdGB , α dCS ), which are the fundamental parameters entering the EdGB and dCS actions. In last two rows of Table I,  Having determined the statistical errors obtained in the single-detector scenario, we can now analyze how these results change when multiple detectors are taken into account. Note that in this case we are assuming that the BBHs are observed simultaneously by all the detectors. This naturally restricts the total number of events to the fraction of binaries detected by LISA which, given the threshold requirement ρ > 15, is approximately 5% of the full population. Therefore, the ensemble of BBHs observed by B-DECIGO, ET, and LV will also be smaller compared to that one analyzed for the stand-alone configuration. Figure 5 shows the joint projected constraints derived using LISA+ET+LIGO/Virgo (empty histograms) and the full network of four detectors (filled histograms), i.e., including B-DECIGO. The median values of the 1σ upper bounds are listed in Table II. For sake of clarity, we also show the constraints for single interferometer, assuming the restricted BBHs ensemble used for the multiband analysis. The latter improves the results achieved by the single detectors, providing stronger bounds. When considering the full network, the errors derived for the dipole correction are very close to those obtained by B-DECIGO alone, which dominates the other interferometers. However, if we neglect the Japanese space detector, the joint LISA-ET-LV analysis improves the bounds on σ B and σ ζ EdGB by more than a factor 20 with respect to the smallest constraints given by LISA. For dCS, the Einstein Telescope plays a more fundamental role for both the network configurations with and without B-DECIGO. Finally, constraints derived for the coupling constants (α EdGB , α dCS ) are also shown in the bottom rows of Table II. The multiband analysis leads to narrower upper limits on the parameters of both quadratic gravity theories, which are factors approximately 10 −2 and 10 −7 smaller than current bounds on EdGB and dCS, respectively [44,47]. These results are also consistent with previous analyses carried out on the evolution of different types of compact binaries in quadratic gravity theories [81][82][83].

VI. CONCLUSIONS
Current and future observations of gravitational wave signals from BBHs carry the full potential of testing General Relativity and tracing modifications of gravity in the relativistic strong-field regime. The constraints on possible non-GR deviations benefit from the synergy of multiple detectors, which greatly improve the accuracy on the measurements of the source's parameters. While LIGO and Virgo are already working as a network of groundbased interferometers, a complementary window will be LISA B-DECIGO ET-D LVC LISA+ET+LV Full network σB 9.0 · 10 −9 1.9 · 10 −11 5.1 · 10 −7 2.7 · 10 −4 3.6 · 10 −10 5.8 · 10 −12 σ ζ EdGB 8.7 · 10 −10 1.9 · 10 −12 4.9 · 10 −8 2.6 · 10 −5 3.5 · 10 −11 5.8 · 10 −13 σ ζ dCS 502 1.7 · 10 −2 2.0 · 10 −1 13 8.0 · 10 −2 3.5 · 10 −3  II. Median values of the distributions for the 68% confidence level bounds on the GR modifications for the multiband analysis. Note that the BBH events analyzed in this case are simultaneously observed by all the interferometers and therefore correspond to a smaller sample than that considered in Table I. The last two rows show the median for the 95% confidence level bounds on the EdGB and dCS couplings.
soon provided from the outer space by LISA. Moreover, third generation detectors like ET, or future satellites such as B-DECIGO, will further enlarge this web of instruments, increasing the accuracy of GW observations. In this paper, we analyze how the different setup of a network composed of LISA, LIGO/Virgo, ET, and B-DECIGO will constrain the fundamental parameters of three scalar-tensor theories, which induce modifications in the waveform at various post-Newtonian orders. Specifically, we consider a generalization of the Brans-Dicke theory and two types of quadratic gravity. We perform a Fisher matrix analysis on the signals produced by a population of stellar-mass, spinning BH binaries. Assuming that GR is the correct theory of gravity, we derive the probability distribution of the upper bounds on the parameters of each alternative theory. Our analysis assesses the effectiveness of multiband gravitational wave observations to put stringent bounds on deviations from Einstein's predictions. Among all detectors, we find that B-DECIGO imposes the strongest constraints on all the theories considered. As expected, LISA and groundbased detectors are able to detect with good accuracy low-and high-frequency modifications of the binary inspiral, respectively.
We then derive joint constraints for the full network of interferometers by combining the results of each instrument. When B-DECIGO is not taken into account, the multiband analysis improves the upper bounds on the GR deviations up to an order of magnitude with respect to the single detector configuration. However, the uncertainties tend to be dominated by the Japanese detector, when the latter is included in the analysis. For a science case, we therefore point out that B-DECIGO plays a fundamental role in testing GR modifications, lead-ing to the tightest bounds on all the alternative theories investigated and therefore supporting the quest for an intermediate-frequency detector between LISA and the ground-based interferometers.
Overall, our results suggest that existing bounds on the three alternative theories of gravity considered in this paper can be improved by several orders of magnitude [19,44,82,84], and that a web of GW detectors able to cover a frequency band within [10 −4 , 10 3 ] Hz, represents a crucial ingredient to achieve precise tests of GR in the strong gravity regime.
Note Added -Recently, a preprint with similar conclusions appeared as an e-print [85].