Scope out multiband gravitational-wave observations of GW190521-like binary black holes with space gravitational wave antenna B-DECIGO

The gravitational wave event, GW190521 is the most massive binary black hole merger observed by ground-based gravitational wave observatories LIGO/Virgo to date. While the observed gravitational-wave signal is mainly in the merger and ringdown phases, the inspiral gravitational-wave signal of GW190521-like binary will be more visible by space-based detectors in the low-frequency band. In addition, the ringdown gravitational-wave signal will be more loud with the next generation (3G) of ground-based detectors in the high-frequency band, displaying a great potential of the multiband gravitational wave observations. In this paper, we explore the scientific potential of multiband observations of GW190521-like binaries with milli-Hz gravitational wave observatory: LISA, deci-Hz observatory: B-DECIGO, and (next generation of) hecto-Hz observatories: aLIGO and ET. In the case of quasicircular evolution, the triple-band observation by LISA, B-DECIGO and ET will provide parameter estimation errors of the masses and spin amplitudes of component black holes at the level of order 1% -- 10%. This would allow consistency tests of general relativity in the strong-field at an unparalleled precision, particularly with the"B-DECIGO + ET"observation. In the case of eccentric evolution, the multiband signal-to-noise ratio by"B-DECIGO + ET"observation would be larger than 100 for a five year observation prior to coalescence, even with high final eccentricities.


Introduction
Among gravitational-wave (GW) events detected by LIGO and Virgo during O1, O2 and O3a runs [1,2], a binary black hole (BBH) merger: GW190521 [3,4] is one of the most striking discoveries. GW190521 is the heaviest BBH merger ever observed, producing the remnant black hole (BH) with the mass of 142 +28 −16 M 1 that can be interpreted as an intermediate mass BH; the source parameters of GW190521 1 We set G = 1 = c with the useful conversion factor 1M = 1.477 km = 4.926 × 10 −6 s. We also assume a "Planck" flat cosmology (when it is needed) with the Hubble constant H 0 = 67.7 km s −1 Mpc −1 , and density parameters Ω M = 0.307 and Ω Λ = 0.694 [5].
arXiv:2101.06402v3 [gr-qc] 27 Dec 2021 (and our notations) are summarized in Table 1. This measurement triggers the intense investigation of GW190521's unique source property. Table 1. Summary of the source parameters of GW190521 as a quasicircular BBH merger, reported by the LIGO-Virgo collaboration [3] based on the BBH waveform model [6]. The symmetric 90% credible interval for each parameter is also quoted. Note that the parameters are written in our notation, and masses are given in the source's rest frame; multiply by (1 + z) to convert to the observer frame. A key element to better understand GW190521 is the precise measurement of the binary parameters (see, for example, Ref. [7] for a possibility of an intermediate mass ratio inspiral). GW190521 is, however, a much heavier BBH system than previously observed GW events, and one of the difficulties here is the short duration and bandwidth of the GW signal that can be observed in the LIGO/Virgo band. In the case of the quasicircular BBH scenario (which is most favored by the LIGO/Virgo analysis [3,4]), the coalescing time and the number of GW cycles at frequency f (in the observer frame) are estimated as t c ∼ 1. 3 (1 + z) − indicating the lack of the GW signals from the sufficiently long inspiral phase. Because of the short duration signal dominated by the merger and ringdown phases, for example, only weak constraints are obtained for the component BH spins and their orientations [4]. Furthermore, even alternative interpretations of the observed GW signal other than massive quasicircular BBH merger in general relativity (GR) would become more relevant; several plausible scenarios are assessed in Section 6 of Ref. [4] by the LIGO-Virgo collaboration. At the same time, the estimation in Eqs.
(1) and (2) suggests a natural way to overcome the hurdle here: observe the inspiral GW signal in the low-frequency band offered by space-based GW detectors. Future GW astronomy in 2030s will utilize the LISA observatory in the milli-Hz band [8] and deci-Hz GW detectors such as B-DECIGO [9]: a prototype GW antenna of the DECIGO mission [10,11] 2 . Figure 1 plots the track of the strain sensitivity curve of GW190521-like non-spinning BBH system, assuming the quasicircular evolution and a simple inspiral-merger-ringdown (IMR) amplitude model given in Eq. (16). 2 Other proposed GW missions in the low-frequency band, including Taiji [12] and TianQin [13] in the milli-Hz band, and MAGIS [14] and TianGO [15] in the deci-Hz band are concisely summarized in, for example, reviews by Ni [16,17].
At the GW frequency f = 0.1 Hz, for example, Eqs. (1) and (2) give the coalescence time t c ∼ 1.1 × 10 5 s and the number of GW cycles N c ∼ 8.1 × 10 3 with the dimensionless characteristic strain h c ∼ 7.6 × 10 −22 (1 + z) 5 where we introduce the source's rest-frame chirp mass M r ≡ ν 3/5 m r t with the symmetric mass ratio ν ≡ m r 1 m r 2 /m r t 2 . In fact, we find in Section 3 that the sky and polarization averaged signal-to-noise ratio (SNR) (whose meaning is momentarily clarified in Section 2) accumulated 5-years before the final coalescence would be ∼ 5.9 × 10 1 in the B-DECIGO band and ∼ 2.7 in the LISA band. These estimations show that the early-inspiral signal of GW190521-like BBHs would be sensitive in the LISA band, and it would be even loud in the deci-Hz band.  Figure 1. Strain sensitivity curves for the ground-based aLIGO and a next-generation (3G) detector: Einstein Telescope (ET), and space-based B-DECIGO and LISA, together with the GW amplitude of a GW190521-like nonprecessing, quasicircular BBH. We use the median values in Table 1 as source parameters. The noise power spectral density (PSD) for each GW detector is given in Section 2.1, and the spectral density of the BBH amplitude is obtained by the nonspinning, IMR amplitude model in Section 2.3. As a reference, we mark with the black dots with 5 years, 1 year, 1 month, 1 day, 1 hour, 1 minute and 1 second before merger time (1), from the left to the right, respectively.
Although the observation with LISA and B-DECIGO alone would provide valuable information of the inspiral GW signal from the GW190521-like BBH system, the true potential of having the low-frequency sensitivity will be revealed only when it is combined with the high-frequency sensitivity in the hecto-Hz band. As illustrated in Figure 1, the late-inspiral and merger-ringdown GW signals of GW190521-like BBH are best detected with aLIGO, Virgo and KAGRA [18]. In addition to ground-based GW observatories that are online, the next generation (3G) of ground-based detectors such as Einstein Telescope (ET) [19,20] (see also Ref. [21] for Cosmic Explorer (CE)) will significantly improve the visibility of GW190521-like BBH systems. We will see in Section 3 that the averaged SNR of the late-inspiral and ringdown signals in the ET band would be ∼ 2.7 × 10 and ∼ 1.5 × 10 2 , respectively. The joint "space + ground" observation across the full GW bands would be therefore the best way to observe the GW190521-like BBH systems; this is the basic idea of multiband GW astronomy.
Soon after the first detection of GW150914 [22], the potential of multiband GW astronomy of BBH systems with LISA and aLIGO was emphasized [23] (see also Ref. [24]). This study was immediately followed up with more detailed analyses. The works include (but not limited to) the improved estimation of source parameters [25,26], tests of GR with high precision [27][28][29][30][31][32], refined event rate estimations [33], probing environment effects [34,35], and new data analysis ideas [36]; all prove the scientific values added by the multiband observation. Among these investigations, Refs. [9,[37][38][39][40][41] demonstrated that the multiband observation of stellar-mass BBH systems will further benefit from having (B-)DECIGO in the deci-Hz band, which naturally bridges the gap between LISA and aLIGO bands.

Goals and organization of this paper
Our purpose in this paper is to explore the prospects for the multiband observation of GW190521-like nonprecessing, quasicircular "intermediate-mass" BBHs. The possibility of multiband observations of intermediate-mass BBHs with deci-Hz GW detectors was pointed out by Amaro-Seoane and Santamaria [42], and Yagi [43]. We consider these observations across the full GW spectrum provided by LISA in the milli-Hz band, B-DECIGO in the deci-Hz band, and aLIGO and ET in the hecto-Hz, looking at two specific aspects of multiband GW astronomy/physics: parameter estimation errors and tests of GR.
We begin in Section 2 by providing a set of our basic tools for the signal analysis in the matched filtering technique: i) the noise PSD in Section 2.1; ii) the GW signal models in Sections 2.2 and 2.3; and iii) the multiband Fisher matrix formalism in Section 2.4. In Section 3, we present the parameter estimation errors of GW190521-like nonprecessing, quasicircular BBH systems using multiband GW observations. They are displayed in Table 2 and 3. Based on the estimated errors, we then examine in Section 4 to what extent the future multiband observations of GW190521-like BBH systems will improve tests of GR; this includes an inspiral-ringdown consistency test (see Table 4), and a simple test to discriminate between the remnant BH based on GR and other remnant compact objects (see Figure 2). We conclude in Section 5 with complications and various effects, which are not covered in this paper due to our assumptions and simplifications.
Appendices contain some additional analysis and information. While a quasicircular BBH merger is the most favored scenario by the LIGO/Virgo analysis [3,4], other alternative scenarios such as an eccentric BBH merger may also be consistent with the observed source for GW190521 [44]. We briefly discuss the multiband visibility of a GW190521-like eccentric BBH in Appendix A. Also, additional noise PSDs of ground/space-based, current/future GW detectors are provided in Appendix B, as a complement of our treatment in Section 2.1.
Throughout this paper, the binary parameters and GW frequencies measured in the source's rest frame are denoted with the index 'r', explicitly distinguished from those in the observer frame.

Method for signal-to-noise ratio and Fisher analysis
In this section, we summarize our methodology of multiband GW data analysis in the specific context of a nonprecessing, spinning, quasicircular GW190521-like BBH system, where the spins are aligned to the orbital angular momentum; an alternative scenario of non-zero orbital eccentricity in GW190521-like BBH system will be discussed in Appendix A. Our simple framework here is largely based on the previous B-DECIGO works [9,38].

aLIGO, ET, B-DECIGO and LISA
For our multiband GW observation, we follow Ref. [38] by considering the four GW observatories: aLIGO and ET in the hecto-Hz band, B-DECIGO in the deci-Hz, and LISA in the milli-Hz band (see also Appendix B for some of other GW detectors in these bands). It should be noted that we shall use the "non sky-averaged" PSD S n ( f ); we will account for the average over the GW detector's antenna pattern function at the level of the waveform.
For aLIGO in the hecto-Hz band, we use Eq. (4.7) in Ref. [45], For a 3rd generation (3G) GW interferometer: ET in the hecto-Hz band, we find Table 1 in Ref. [46] as For B-DECIGO in the deci-Hz band, we use Eq. (20) in Ref. [38], originally proposed by Nakamura et al. [9]: For LISA in the milli-Hz band, we use Eq. (1) of in Ref. [47], which is based on the 2018 LISA Phase-0 reference design parameters. It reads 3 where S LISA c ( f ) represents the effective PSD due to the unresolved galactic binaries; the explicit expression is given in Eq. (14) of Ref. [47] and we have assumed a 4-year mission for Figure 1. 3 Our expression S LISA n (7) is smaller than Eq. (1) in Ref. [47] by a overall factor of 5; the difference arises simply because Eq. (7) does not account for the sky-averaging.

Waveform models
We employ as our BBH waveform model the frequency-domain, "restricted" waveform in the stationary phase approximation for the inspiral phase, and the frequency-domain, single-mode waveform for the ringdown phase. We shall restrict our waveform model to these two phases to simplify our analysis as far as possible; the complete IMR treatment at the level of waveform (using, for example, "effective-one-body" (EOB) approach [48,49] and "phenomenological" (IMRPhenom) model [50,51]) will be left for future work.
The GW waveform from a BBH inspiral in the frequency domain has the well-known form (see, for example, Ref. [52])h where A is the "Newtonian" amplitude averaged over all sky positions and binary orientations (see, for example, Ref. [53]), so that A ≡ 2 5 The waveform's frequency-domain phase Ψ Insp ( f ) in the post-Newtonian (PN) approximation is given by (see, for example, Ref. [54]) where v ≡ (πm t f ) 1/3 is the PN parameter (in terms of the observer-frame total mass), and t c and Ψ c are the time and phase at coalescence. The phase terms ∆Ψ are the 3.5PN spin-independent, point-particle contributions derived in Ref. [55] and the 3.5PN spin-dependent, point-particle contributions that include linear spin-orbit [56,57], quadratic-in-spin [58] and cubic-in-spin [59] effects, respectively. The remaining phase term ∆Ψ BH−tidal 3.5PN is related to the tidal response of a spinning BH as a finite-size body, i.e., BH-absorption corrections such as the GW energy and angular momentum fluxes down to the horizons and the associated evolution of the BH itself [60][61][62].
Meanwhile, limited to only fundamental (n = 0), = 2 = m mode, the time-domain, single-mode ringdown waveform measured at a GW observatory is written as [63][64][65] h where t 0 and φ 0 are the initial time and phase of the ringdown, respectively, and we have ignored the overall amplitude so that Eq. (11) is not normalized; the initial ringdown amplitude may be determined by matching the ringdown GW waveform to the merger one (refer to, for example, Eq. (16) below). When the final remnant object is a Kerr BH, the central frequency f c and the quality factor Q are given in terms of quasinormal-mode (QNM) frequencies ( It should be noted that f QNM = f r QNM (1 + z) −1 here is given in the observer frame. The Fourier transforms of Eq. (11) provide the corresponding frequency-domain waveform, which takes the form (here we follow the convention of Ref. [63]) We note that this frequency-domain GW waveform is not normalizedeither 4 , following the approach presented in Appendix B of Ref. [64]. The ringdown waveforms (11) and (13) require input data for the QNM frequency of the final remnant BH. Because we do not consider the complete IMR phase at the level of waveforms, we employ the numerical-relativity (NR) remnant fitting formulas (see, for example, Refs. [68][69][70][71][72][73] and references therein), from which the final mass M f and spin S f of the remnant BH are consistently inferred for a given initial BH masses m 1,2 and dimensionless spin parameters χ 1,2 (in the inspiral phase) as 5 We then generate the accurate numerical data of f QNM of the inferred remnant BH with the Black Hole Perturbation Club (B.H.P.C.) code [75] to obtain f c and Q 6 . In practice, it is also convenient to present f c and Q by means of a compact analytical formula. Such a formula for the ( = 2, m = 2, n = 0) mode is obtained in Ref. [65] (by performing fits to the numerical QNM frequency data), and it reads with f 1 = 1.5251, f 2 = −1.1568, f 3 = 0.1292, q 1 = 0.7000, q 2 = 1.4187 and q 3 = −0.4990.

Signal-to-noise ratio
We estimate the SNR of a complete IMR GW signal, making use of the simple frequency-domain, IMR "amplitude" model [9]. This model is motivated by the waveform amplitude in the IMRPhenomB model [80] and it is given by 4 The maximized SNR over the initial ringdown phase φ 0 has been discussed in Ref. [67]. 5 The mass and spin of the remnant BH in the remnant formulas are derived in the isolated horizon framework (see, for example, Ref. [74]), not obtained from the ringdown GWs. The latter is used only for checking the internal consistency of the formulas. 6 For accurate numerical data of f QNM , see also Ref. [76], Emanuele Berti's "Ringdown" website [77], Ref. [78] and the Black Hole Perturbation Toolkit [79].
where the overall constant A is chosen to be the (averaged) GW signal's amplitude in the inspiral phase (9). We set f max = 1/[6 3/2 π(1 + z)m r t ] as the GW frequency at the innermost stable circular orbit (ISCO) of a test particle in the Schwarzschild spacetime with the total mass of BBH m r t , and f r R and f r I are the real and imaginary parts of the QNM frequency (i.e., f r QNM = f r R + i f r I ) of the fundamental (n = 0), = 2 = m mode, which are determined by the final mass and spin of the remnant BH; the redshift dependence appears because the model parameters (m r t , f r R , f r I ) are all given in the source's rest frame (while the GW frequencies f and f max are in the observer frame). The averaged SNR (in the above sense) can be then obtained by 7 where S n ( f ) is the noise PSD, and the frequency range [ f in , f end ] is determined by the GW detector with which we observe the GW signals.
Note that the amplitude spectral density (i.e., the square root of the PSD of the source amplitude) is | in this model, and its track for the GW190521-like BBH is plotted in Figure 1 with the (square root of) the noise PSD S n ( f ), assuming the median values of Table 1 and the remnant formulas of Ref. [68].

Multiband Fisher analysis
We approximate the variance (i.e., uncertainty squared) associated with the measurement of a set of signal parameters, making use of the standard Fisher matrix formalism. The Fisher information matrix for a single-band GW detector is defined by whereh( f , θ) is the frequency-domain GW signal described by the set of parameters θ, and θ 0 are the best-fit values of the binary parameters. The bracket defines the noise-weighted inner product over the with asterisk ' * ', denoting the complex conjugation. The inverse Fisher matrix defines the corresponding variance-covariance matrix Σ ab ≡ (Γ ab ) −1 . In the limit of suitably high SNR [83], the variance of the parameter θ a is given by 7 For the given Fourier transform of a GW signalh( f ), the SNR can be written in terms of eitherh( f ) itself with the dimension of 1/Hz, the spectral density of the source amplitude [81]: In the case of the multiband analysis (to combine the information from, for example, aLIGO + B-DECIGO), we simply construct a multiband SNR and Fisher matrix by adding individual (averaged) SNRs and Fisher information matrices for each GW detector: where ρ I ave and Γ I ab are the averaged SNR and the Fisher matrix for the I-th detector. The multiband variance-covariance matrix is defined by and the variance of θ a is then obtained by σ 2 a = Σ aa tot .

Parameter estimation errors via multiband observation
In this section, we summarize the parameter estimation errors of the inspiral and ringdown phases for a nonprecessing, spinning, GW190521-like BBH system, using the multiband GW network (LISA, B-DECIGO, aLIGO and ET) detailed in Section 2.1. We follow Refs. [38,66] in our treatment of the Fisher matrix calculation for BBH GW signals 8 , and we continue to neglect the contribution from the merger GW signal 9 ; for the reason that one cannot separate the inspiral, merger and ringdown phases cleanly in the strict sense although we have presented the simple IMR amplitude model in Eq. (16). When the merger contribution is introduced into the inspiral or ringdown signal analysis, it causes some bias in the inspiral or ringdown parameter estimation. For example, the merger-ringdown waveform is parametrized by the binary parameters, not solely by the remnant BH parameters after the merger. Similarly, in the inspiral-merger waveform, the "late" merger phase can be described by the overtones of QNMs [85], i.e., the remnant BH parameters. Thus, we shall perform the inspiral-ringdown consistency test of GR without the merger contribution in Section 4.1.

Setup of Fisher analysis
We set the default frequency interval of each GW detector [ f low , f up ] as [10.0, 3.0 × 10 3 ] Hz (aLIGO), [2.0, 3.0 × 10 3 ] Hz (ET), [0.01, 1.0 × 10 2 ] Hz (B-DECIGO), and [1.0 × 10 −4 , 1.0] Hz (LISA), respectively. We shall adopt the T obs = 5-year observation time, and assume that the binary merges at the GW frequency of the Schwarzschild ISCO, f ISCO = 1/[6 3/2 π(1 + z)m r t ]. In this setup, the minimum frequency of the GW signal is where M r is the chirp mass in the source's rest frame, normalized to that of GW190521. Therefore, the GW signal observed by each GW detector is truncated at the corresponding initial frequency f in ≡ max( f min , f low ) as well as the end frequency f end ≡ min( f ISCO , f up ); recall Figure 1. 8 The setup here is slightly different from Ref. [38]; i) we will assume the 5 yr observation, rather than 4 yr observation, and ii) we will use the new LISA sensitivity curve proposed by Ref. [47] and displayed in Eq. (7), not the earlier eLISA sensitivity curve presented in Ref. [84].
The parameters of the inspiral waveform (8) are where we define the symmetric and anti-symmetric combinations of BH spins by χ s ≡ (χ 1 + χ 2 )/2 and χ a ≡ (χ 1 − χ 2 )/2 with the component (aligned) BH spins χ 1,2 ≡ | χ 1,2 | ≡ |S 1,2 |/m 2 1,2 . At the same time, the parameters of the ringdown waveform (13) are It should be noted that the amplitude parameters are left out from the set of our independent parameters both in Eqs. (24) and (25). They are entirely uncorrelated with other parameters θ a because the variance-covariance matrix Σ a b gives the variance σ 2 ln A = ρ −2 ave and the correlation c ln A, a ≡ Σ ln A, a /(σ ln A σ a ) = 0 for the inspiral GW signal (8) (see, for example, Ref. [86]), and similar for the ringdown GW signal (if we explicitly introduce the amplitude to the normalized waveform of Eq. (13)). For simplicity, we consider that all (other) parameters are unconstrained.
The best fit values of inspiral parameters are given by the median values in Table 1 with t c = 0.0 = Ψ c , while those of the ringdown parameters are assumed to be φ 0 = 0.0 = t 0 10 , and Note that these values are different from the mass and spin (M f , χ f ) = (142 +28 −16 M , 0.72 +0.09 −0.12 ) of the remnant BH reported in the LIGO/Virgo GW190521 detection paper [3]; we have assumed that the individual spins of GW190521-like BBH are nonprecessing, and completely aligned to the orbital angular momentum for simplicity while the observed GW190521 is actually considered to be a precessing BBH. The associated QNM frequency data of this remnant BH is then generated by the B.H.P.C. code [75], yielding the results in Eq. (26).
The root-mean-square of parameter estimation errors scales like the inverse of SNR, ∼ 1/ρ, and it depends on both ρ and bandwidth over which the SNR is accumulated. To see the benefit of having wider bandwidth due to the multiband observation, finally, we introduce the normalized root-mean-square error as and shall display our error estimations in terms of δθ with the total averaged SNR ρ tot ave accumulated over multi-frequency bands. 10 The parameter estimation errors are independent of the value of the initial time t 0 . In our frequency-domain, single-mode ringdown waveform in Eq. (13), the t 0 dependence is factorized as e 2 iπ f t 0 and it does not contribute to the noise-weighted inner product in Eq. (19). On the other hand, the estimation errors depend on the initial phase φ 0 weakly [64].

Result: Inspiral phase
In Table 2, we present the parameter estimation errors of mass parameters (m, ν) and spin parameters (χ s , χ a ) for the GW190521-like BBH inspiral in various combinations of ground/space-based GW observatories (but suppressing those of t c and Ψ c ). Here, it should be noted that only the inspiral phase is analyzed here, and the merger-ringdown phase which contributes to the SNR for ground-based observatories, is ignored ; recall the footnote 9. Therefore, the SNR for aLIGO quoted here is much smaller than the observed LIGO/Virgo network SNR of 14.5 [3].
In the single-band case, as seen in Figure 1, the inspiral GW signal of the GW190521-like BBH is best observed by B-DECIGO because it can cover both the early (1 yr before merger) and late (around the ISCO) phases. However, B-DECIGO observation alone is not enough benefit to discern the spin parameters.
In the multiband cases, thanks to the wider bandwidth, the observation with B-DECIGO and ET gives a factor of 2 improvement in all the parameter estimation, even if we observe only the inspiral GW signal. This is further refined if we combine the data from LISA, forming a triple-band network (although the estimated LISA SNR of ∼ 2.68 is likely too small compared to the 'realistic' detection SNR threshold of ∼ 15 [87]). Importantly, the normalized errors of BH spins (δχ s /χ s , δχ a /χ a ) = (2.52, 5.71 × 10 2 ) in this best case imply that the magnitudes of individual BH spins σ χ can be recovered with the fractional statistical errors (i.e., now accounting for the total multiband SNR of ρ tot ave = 6.52 × 10 1 to the uncertainties; recall Eq. (28)) (σ χ 1 /χ 1 , σ χ 2 /χ 2 ) = (2.43 × 10 −1 , 2.51 × 10 −1 ). The joint LISA, B-DECIGO and ET observatories would be only viable network to measure all the binary parameters of GW190521-like BBH system in the inspiral phase, including BH's component spins. Table 2. Parameter estimation errors of mass parameters (m, ν) and spin parameters (χ s , χ a ) for the GW190521-like quasicircular BBH inspiral, normalized to the total multiband SNR ρ tot ; for example, the result of "BD" in "BD + aLIGO" is normalized to ρ tot = 5.93 × 10 1 . We assume that the true binary parameters are given by the median values of Table 1, i.e., m t = 151.0 M , ν = 0.246, χ s = 0.71 and χ a = 0.02. Since the source location can be determined well by B-DECIGO (BD here) [9,39], we fix 1 + z = 1.82. Note that only the inspiral SNR of aLIGO is 1.76, which is too small to give meaningful estimation errors. Also any normalized estimation errors δθ > 1.0 × 10 6 are discarded from this

Result: Ringdown phase
In Table 3, we show the parameter estimation errors of the central frequency f c and quality factor Q (but suppressing those of t 0 and φ 0 ) for the ringdown phase of the GW190521-like BBH in ground-based observatories. For simplicity, we estimate the ringdown amplitude and associated SNR via the IMR amplitude model in Eq. (16) 11 . We see that ET will be able to measure the QNM frequency of the remnant BH with the statistical error ∼ 10 −3 . At the same time, however, the difference in the normalized errors δf c and δQ between aLIGO and ET observations is not so evident. We speculate that it arises from the difference in the spectrum shapes (not the overall amplitudes) because the ringdown waveform (13) is narrow banded in the frequency domain.

The implications for tests of GR via multiband observation
In this section, we explore what extent the multiband observation of the GW190521-like BBH system discussed in Section 3 could improve tests of GR. A handful of tests have been already formulated and performed with merging BBH systems (see, for example, Refs. [88][89][90][91]), and we follow the (very) simple tests proposed by Nakano et al. [66,92].

A consistency test of GR with the inspiral and ringdown GW signals
One possible test of GR with a BBH system is to establish the consistency of the mass and spin of the final remnant BH determined by two different parts of the GW signals. Thanks to the recent advancement in NR simulations of BBH systems [93][94][95] (see also Refs. [96][97][98][99][100]), one can infer these values from the initial component masses and spins measured from the inspiral GW signal (in the low-frequency band), making use of the NR fitting formulas for the remnant properties of the final BH; recall Section 2.2. At the same time, they are directly estimated from the succeeding merger-ringdown GW signal (in the high-frequency band). This type of test is now known as the "IMR consistency test" [101,102]. By formulation, multiband observations of heavy BBH mergers such as GW190521-like BBH systems will be "golden binaries" [103] of such an IMR consistency test.
Given that the early inspiral and late ringdown GW signals will be best observed in a different frequency band (such as "B-DECIGO + aLIGO" network etc.), we here perform the multiband version of the "inspiral-ringdown consistency test" formulated by Nakano et al. [66] (see also Refs. [103,104]), solely using the inspiral and ringdown parts, and test the consistency of GR across the merger part, which is a highly dynamical phase in a strong-field regime. We estimate the statistical errors on M f and χ f from the inspiral GW signal by using the (normalized) statistical errors δθ in Table 2 and applying a standard variance propagation of non-linear functions to the specific NR remnant formulas ("UIB formulas") [71] publicly available in LALInference [105,106] 12 . The errors from the ringdown GW signal are estimated from Table 3, through the dependence of f c and Q on (M f , χ f ) 13 . 11 We should note that this approximation is likely too raw because the ringdown amplitude strongly depends on the starting time t 0 of the ringdown phase, which is difficult to determine in practice. 12 We ignore the systematic bias due to our specific choice of the NR remnant formulas, for simplicity. See, for example, Refs. [99,102] for details. 13  In Table 4, we present the parameter estimation errors of the remnant mass M f and spin χ f from both the inspiral and ringdown phases. For the reference, we also present the same errors of the GW150914-like BBH with the redshifted component masses (m 1 , m 2 ) = (30 M , 40 M ), spin magnitudes (χ 1 , χ 2 ) = (0.9, 0.7), and the luminosity distance D L = 0.4 Gpc (i.e., the redshift z ∼ 0.085) 14 . We see that the "LISA + B-DECIGO + ET" observation allows us the IMR consistency test at the sub-percent precision both for the GW150914-like, and GW190521-like BBH systems. The result of the GW190521-like BBH system is slightly worse than the case of the GW150914-like BBH system because the luminosity distance of GW190521 (D L = 5.3 Gpc) is much larger than that of GW150914 (D L = 0.4 Gpc) [22,107,108]. Also, the GW190521-like BBH system has the central frequency f c ∼ 85 Hz lower than that of the GW150914-like BBH system ( f c ∼ 350 Hz), missing the most sensitive frequency of aLIGO and ET ∼ 250 Hz. If a GW190521-like BBH was observed at the same distance as GW150914, the accuracy of its inspiral-ringdown consistency test would be at the level of ∼ O(0.01%). Table 4. Parameter estimation errors of the remnant mass and spin of GW190521-like and GW150914-like BBH systems, using only the ringdown GW signal as well as the inspiral GW signal inferred via the NR remnant fitting formulas; the best fit values (M f /m, χ f ) ∼ (0.904, 0.883) for the GW190521-like BBH (recall Eq. (27)) and (M f /m, χ f ) ∼ (0.891, 0.898) for the GW150914-like BBH. Here, BD is an abbreviation for B-DECIGO. The results in the single band are normalized to the SNR for a given GW observatory. In the multiband case they are normalized to the total SNR of ρ LISA+BD+ET .

A simple test of the remnant compact object with quasinormal modes
Another simple test of GR is to bracket whether the remnant object should be a BH predicted by GR or not, making use of the parameter estimation errors of QNM frequencies ( f QNM = f R + i f I ) obtained from ringdown GW signals [66]. Figure 2 plots the 1σ, 2σ, 3σ, 4σ, and 5σ error contours on the fundamental (n = 0), = 2 = m mode of the QNM frequency in the ( f R , f I ) plane, in the case of the GW190521-like BBH system observed by aLIGO (left) and ET (right). The errors are estimated through the results in Table 3 about the parameter estimation error on the ringdown GW signal (after t 0 and φ 0 being marginalized out), and the outermost contour in each panel shows the 5σ error. This black lines in each panel depict the Schwarzschild limit of | f I |/ f R , which may be obtained by setting χ f = 0 in Eq. (15) with Eq. (12) 15 , The key point of this test is that the top-left side of the black line becomes the prohibited region in GR, i.e., for the (n = 0), = 2 = m mode, the QNM frequencies ( f R , f I ) of any rotating Kerr BHs must sit in the bottom-right side of the black line.
In the aLIGO case, due to the low SNR (= 10.8), the parameter estimation errors already go beyond the Schwarzschild limit at the 3σ level, and we cannot confirm whether the remnant object is a BH predicted by GR with the 5σ level. For such low SNR events, the "coherent mode stacking method" [109] will be useful. On the other hand, thanks to the high SNR (= 147) in the ET case, this simple test can necessarily confirm that the remnant object is a GR-predicted BH 16 . Figure 2. Simple remnant test in the ( f R , f I ) plane. The QNM frequency is given in the observer frame, and we assume the GW190521-like ringdown GW signals with (M f /m, χ f ) ∼ (0.904, 0.883) (27). The figure shows the parameter estimation error with the aLIGO noise curve (left, SNR = 1.08 × 10 1 ) and the ET one (right, SNR = 1.47 × 10 2 ). The black lines are the Schwarzschild limit of | f I |/ f R ≈ 0.236. The colored ellipses show the contours of 1σ, 2σ, 3σ, 4σ, and 5σ. To obtain these two dimensional plots, the time and phase parameters (t 0 , φ 0 ) have been marginalized out.

Summary and Discussion
This work underlines the multiband observation of the GW190521-like nonprecessing, quasicircular "intermediate-mass" BBH with LISA and B-DECIGO in the low-frequency band, combined it with aLIGO and ET in the high-frequency band. Our first result of the parameter estimation errors was displayed in Tables 2 and 3; the statistical errors of the binary's mass parameters by B-DECIGO observation will be ∼ 10 −2 , and the multiband observation with LISA, B-DECIGO and ET will further improve them to a factor of 2, even allowing the statistical errors of component BH spins at the accuracy level of ∼ 10 −1 . Based on the ringdown analysis, ET will measure the QNM frequency about O(0.1)% precision. Our second result of GR tests was presented in Table 4 as well as Figure 2. We showed that the multiband observation of the GW190521-like BBH system by LISA, B-DECIGO and ET would perform the inspiral-ringdown consistency test at the level of percent-precision.
The main point of our analysis is that there is the principal advantage of measuring GW190521-like BBH systems using the full GW spectrum from milli-Hz to deci-Hz bands, and to hecto-Hz band. We expect our findings will motivate further investigation on prospects for multiband observation of the GW190521-like BBH systems, as much alike the prototypical GW150914-like BBH systems [22].
Nevertheless, we emphasize that our study was performed with (very) simple methods. Therefore, our results should be only indicative and tentative. In the remainder of this section, we discuss what works remain to be done to refine our analysis in future, so that we can eventually make a strong scientific cases of multiband GW astronomy/physics.

Assessment of prospects
First, our methodology in Section 2 should be replaced with more modern, sophisticated approaches to the GW data analysis; it will include (but not limited to), for example, the use of complete IMR waveforms such as EOB approach as well as IMRPhenom models (see, for examples, Ref. [110] and references therein), full-fledged Bayesian posterior based techniques (see, for examples, Refs. [105,111,112]), and more 'realistic' noise and waveform models that account for the sky-location of sources as well as orbital configurations of B-DECIGO [11] and LISA [8].
Second, our target GW190521-like BBH system was restricted to the nonprecessing, spinning, quasicircular configuration; although there is large uncertainty, the spin precession of GW190521 is estimated as nonzero [3,4]. Future studies are therefore needed to concern both the spin precession and orbital eccentricity of the BBH system: we will briefly discuss the eccentric, non-spinning BBH system in Appendix A. In general, the BH's spins in the precessing binary have not only their magnitudes but also orientations, which can be described by, for examples, the effective inspiral spin parameter χ eff (related to the components aligned with the orbital angular momentum), and the precession spin parameter χ p (related to the components in the orbital plane for the inspiral GW waveform). Adding the orbital eccentricity to the (precessing) BBH system will be fully generic, and the waveform modeling becomes (much) more complicated. Despite that challenge, there is a considerable development of the analysis on fully generic binary systems [113].
Related with the point mentioned above, it should be noted that we have ignored subdominant ( = 2, |m| = 2) harmonics in our GW waveform model. They are more notable in the observed signal when system's mass ratio becomes smaller (such as the analysis of GW190814 [114] with the mass ratio q = 0.112 +0.008 −0.009 ). With the subdominant harmonics, one can access the source orientation to reduce uncertainty in the distance estimation [115].
Third, we should note that there are two main hurdles to analyze the ringdown GW signals; the low SNR with aLIGO, and the starting time of the ringdown phase that is a priori unknown to the whole observed GW signal (see, for examples, Refs. [116,117] for discussions on the starting time). Our simple analysis is carried out with the single-mode waveform model (11) as the template to analyse the ringdown phase, assuming that the starting time of the ringdown phase to be t = t 0 which corresponds to the f = f R , i.e., just after the end of the merger phase. While this choice does not affect our results of parameter estimation errors, in practice, the best fit values (that should be obtained rather than assumed in the context of the full parameter estimation against the raw GW data) can be biased if one assumes the earlier starting time in the analysis (see, for example, Figure 5 in Ref. [88]). Although one may delay the starting time to avoid the bias in the parameter estimation, the SNR becomes much lower than the expected SNR with the assumption of t = t 0 .
These obstacles will be overcome if one includes higher overtones (n > 0) into the ringdown GW analysis [85], for which a much larger SNR than the single-mode analysis will be expected. Indeed, a superposition of overtones (n > 0) in a single harmonic mode will be observed [85,[118][119][120] in the high SNR events (with ET). Another step-functional improvement of the ringdown GW analysis will be offered by using completely different signal-analysis methods than the traditional matched filtering analysis, which may not be always optimal for the ringdown GW signal. There is ongoing work to assess the improvement due to new techniques for the ringdown GW signal analysis (such as Hilbert-Huang transformation, autoregressive modeling, and neural networks) [121,122].
Fourth, there are other GR tests that were not covered in Section 4, but can be greatly improved using the multiband observation of GW195021-like BBH systems. Reference [91] performs various tests of GR with the BBH events in GWTC-2 (see also the reviews by Carson and Yagi [123] about the current and future test with GWs), including i) the IMR consistency test between the inspiral and postinspiral phases divided at some cutoff frequency; ii) constraining deviations from GR with parametric deformations to a predicted GR waveform model, iii) "BH spectroscopy" [124,125] with ringdown GWs which contains two (or possibly more) QNMs [126,127], and so on.
The test i) is similar to our inspiral-ringdownconsistency test directly using information of the merger phase, and makes the most of the GW waveform. It has been pointed out that the values of this test with GW150914-like BBH systems will be maximized in the multiband observation [28,29,128]. Our result suggests that the same will be true for GW190521-like "intermediate-mass" BBH systems, too. Similarly, Gupta et al. [31,32] have showed that the multiband observations of stellar-and intermediate-mass BBHs with LISA and 3G detectors will be only workable way to carry out the most general version of test ii). Adding B-DECIGO (or any other planned GW detectors) in the deci-Hz band to the multiband analysis, we expect the precision of this test will be unprecedented.
Based on the generalized likelihood ratio test [129], the test iii) is performed (see, for examples, Section 9.5 in Ref. [130] and references therein, and also Ref. [131] for the future O5 era). This test with aLIGO and Advanced Virgo alone is quite challenging to have any conclusive result, simply because of low SNRs and larger parameter estimation errors in the ringdown phase. In the 3G era, two (or possibly more) QNMs will be measurable [126,127]. It is also helpful to use the multiband observation in order to optimize ground-based detectors via the forewarnings from the low-frequency, LISA band [132]. In either cases, the multiband observation with B-DECIGO and ET will give them additional advantage of being able to perform the best test of Kerr hypothesis of remnant BHs via BH spectroscopy.

Appendix A. Signal-to-noise ratio of GW190521-like eccentric BBH systems
Throughout the bulk of this paper, we have looked at the GW190521-like BBH system under the assumption of a quasicircular BBH merger. While the quasicircular evolution of GW190521 is totally consistent with the LIGO/Virgo observation [3,4], due to the lack of the inspiral GW signal long enough, it appears that alternative scenarios, for example, GW190521 as an eccentric BBH merger also become relevant. Indeed, Gayathri et al. [44] demonstrated that observed GW190521 data could be explained as an equal-mass, highly-eccentric (e = 0.7) BBH system 17 , and the estimated source parameters are quite different from those derived from the quasicircular BBH scenario (recall Table 1 (see Ref. [44] for the spin parameters). The possibility of GW190521 with nonvanishing eccentricity is also pointed out by Refs. [134,135]. Like the quasicircular case, the multiband observation of eccentric GW190521-like BBH systems will once again help in distinguishing these two scenarios. Assuming a quadrupole GW generation from a Newtonian Kepler orbit [136,137], the typical coalescing time t ecc c and the characteristic strain amplitude h ecc c, n of the n-th harmonic are (see, for example, Section 4.1 of Maggiore's text [138] as well as Ref. [139]) where t circ c and h circ c are corresponding circular-inspiral results given in Eqs.
(1) and (3), respectively, and the function g(n, e) will be defined momentarily. We see that the early inspiral phase of the eccentric binary is well within the B-DECIGO and LISA bands, too. Specifically, Holgado et al. [140] (see also the works by Amaro-Saoane [42,139,141]) pointed out that having deci-Hz GW observatories such as B-DECIGO, MAGIS [14] and TianGO [15] will be a key element to observe the eccentric inspiral GW signals in multiband networks because the GW signals may skip the LISA band entirely; h ecc c, n can be suppressed by a function of g(n, e) significantly, depending on its eccentricity in the LISA band.
To better understand the visibility of eccentric BBH systems with aLIGO, ET, B-DECIGO and LISA, let us estimate the SNR of non-spinning, eccentric BBH inspirals accumulated in each band; see also Ref. [142] for a recent review about waveform families for the eccentric binary systems. We also see various active works on eccentric waveform approximants [143][144][145][146]. The squared SNR averaged over the all-sky positions and binary orientations may be written as [23,53,81,147] 18 where S( f ) is the noise PSD for a given GW detector, and is the frequency of the harmonic in the observer frame, defined by the source's rest-frame frequency f r orb of the Kepler orbit with the redshift z. The dimensionless characteristic strain h c, n of the n-th harmonic is [147] h c, n ≡ where D L is the luminosity distance. Assuming the quadrupole formula applied to the Kepler orbit (with the chirp mass M in the observer frame) [136,137], is the emitted GW energy per unit frequency f r n at the n-th harmonic measured in the source's rest frame, and we define g(n, e) ≡ n 4 32 with the Bessel functions of the first kind J n (x) (n: integer); see, for example, Ref. [148] for the derivation of Eq. (A5).
For an inspiraling eccentric BBH, the evaluation of SNR through Eq. (A2) requires the knowledge of slowly-evolving orbital eccentricity e and the frequency f n in time, under the gravitational radiation losses. Again, in the quadrupole formalism, it is given by [136,137] with reference eccentricity e 0 and orbital frequency f orb, 0 (in the observer frame). One can set these constants by the values at the last stable orbit (LSO) of the eccentric geodesic (of a test particle) in the Schwarzschild geometry [149]: namely, e 0 = e LSO and with the total mass m r t in the source's rest frame. The frequency evolution (A8) is therefore completely determined with a given single parameter e LSO .
We compute the averaged SNR ρ ave given in Eq. (A2) for the GW190521-like, non-spinning, eccentric BBH with the source-frame masses (m r 1 , m r 2 ) = (102 M , 102 M ) and the luminosity distance D L = 1.9 Gpc (i.e., the redshift z ∼ 0.35), which mimics the eccentric BBH merger obtained by the NR simulations in Ref. [44]. We assume the five year observation prior to the final merger determined by Eqs. (A1) and (A9), and apply the same setup described in Section 3.1 to each n-th harmonic of GW strains. Because the frequency evolution (A8) is expressed in term of eccentricity, in practice, we change the integration variable of Eq. (A2) from f n to e for the computational efficiency, making use of d f n = n|d f orb /de| de with (see, for example, Ref. [150]) f orb e F(e) (1 − e 2 ) 9/2 1 + 121 304 e 2 . (A10) Also, the infinite summation over the harmonics n in Eq. (A2) is truncated at some finite value of n max ; we chose n max = O(10 4 ) so that the resultant SNRs are guaranteed to have at least 3 significant digits 19 .
In Table A1, we summarize the averaged SNRs ρ ave of the GW190521-like, non-spinning, eccentric BBH accumulated in each GW band, for sample values of final eccentricities at the last stable orbit e LSO = {10 −6 , 10 −3 , 0.1, 0.4, 0.6}. For references, the 'initial' eccentricities at f in for each band are also listed. It is approximated by solving Eq. (A8) for e at the detector's initial GW frequency of the second harmonics f orb = f in /2. We found three main results: i) B-DECIGO and ET always accumulate a total SNR greater than at least 10, independent of the value of the final eccentricity e LSO ; ii) LISA has the detectable SNR only when e LSO < 10 −3 . That is, the inspiral GW signal from the GW190521-like BBH that has the high-eccentricity in the aLIGO band would entirely skip the LISA band; iii) the SNR with B-DECIGO becomes bigger when e LSO becomes smaller, while the SNR with ET and aLIGO shows the opposite behavior. Therefore, the multiband observation could provide much louder SNR than the single-band SNR across the full range of e LSO . Using Eq. (21) with the result in Table A1, we find that the multiband SNR with B-DECIGO and ET are ∼ 180 for e LSO = {10 −6 , 10 −3 }, ∼ 140 for e LSO = {0.4, 0.6} and ∼ 100 for e LSO = 0.1. Table A1. Averaged SNRs ρ ave of the GW190521-like, non-spinning, eccentric BBHs accumulated in each band five years prior to coalescence, for given values of final eccentricity e = e LSO . We assume the source-frame component masses (m r 1 , m r 2 ) = (102 M , 102 M ), and the luminosity distance D L = 1.9 Gpc (i.e., the redshift z ∼ 0.35) [44]. The values in parentheses indicate the 'initial' eccentricity estimated from the initial GW frequency ( f in defined in Section 3) of the second harmonic mode in each detectors. Note that aLIGO's 'initial' eccentricity when e LSO = 0.6 is not displayed because the second harmonic mode is not detectable in this case.
SNR and eccentricity at f in e LSO = 10 −6 e LSO = 10 −3 e LSO = 0. Finally, although the parameter estimation is not the main focus here, we briefly discuss the potential accuracy of the eccentricity measurement. In the small-eccentricity and high-SNR limit, the orbital eccentricity may be measured within the fractional error [153] (see also Ref. [154]) from Eq. (2) in Ref. [153] by using a rough approximation with the quasicircular amplitude and the eccentric phase, where δê 0 denotes the parameter estimation error of e 0 (at the GW frequency f 0 for n = 2) normalized by the SNR. Here, the power α is due to the approximation of the noise PSD by the power law, S n ∼ f 2α (assuming α > −2/3), and B-DECIGO may have α = 1, for example. This estimator implies that B-DECIGO would be able to precisely measure the eccentricity of GW190521-like BBH systems. Because we find that the 'B-DECIGO + ET' combination always provides the multiband SNRs larger than ∼ 100 independent of the value of e LSO , one might expect that this combination best observes GW190521-like BBHs over the full range of eccentricity, helping to understand the population properties of BBH mergers [155]. We will explore this possibility in future work.