Multiband Gravitational-Wave Astronomy: Observing binary inspirals with a decihertz detector, B-DECIGO

An evolving Japanese gravitational-wave (GW) mission in the deci-Hz band: B-DECIGO (DECihertz laser Interferometer Gravitational wave Observatory) will enable us to detect GW150914-like binary black holes, GW170817-like binary neutron stars, and intermediate-mass binary black holes out to cosmological distances. The B-DECIGO band slots in between the aLIGO-Virgo-KAGRA-IndIGO (hecto-Hz) and LISA (milli-Hz) bands for broader bandwidth; the sources described emit GWs for weeks to years across the multiband to accumulate high signal-to-noise ratios. This suggests the possibility that joint detection would greatly improve the parameter estimation of the binaries. We examine B-DECIGO's ability to measure binary parameters and assess to what extent multiband analysis could improve such measurement. Using non-precessing post-Newtonian waveforms with the Fisher matrix approach, we find for systems like GW150914 and GW170817 that B-DECIGO can measure the mass ratio to within $<0.1\%$, the individual black-hole spins to within $<10\%$, and the coalescence time to within $<5\,$s about a week before alerting aLIGO and electromagnetic facilities. Prior information from B-DECIGO for aLIGO can further reduce the uncertainty in the measurement of, e.g., certain neutron star tidally-induced deformations by factor of $\sim 6$, and potentially determine the spin-induced neutron star quadrupole moment. Joint LISA and B-DECIGO measurement will also be able to recover the masses and spins of intermediate-mass binary black holes at percent-level precision. However, there will be a large systematic bias in these results due to post-Newtonian approximation of exact GW signals.

An evolving Japanese gravitational-wave (GW) mission in the deci-Hz band: B-DECIGO (DECihertz laser Interferometer Gravitational wave Observatory) will enable us to detect GW150914-like binary black holes, GW170817-like binary neutron stars, and intermediate-mass binary black holes out to cosmological distances. The B-DECIGO band slots in between the aLIGO-Virgo-KAGRA-IndIGO (hecto-Hz) and LISA (milli-Hz) bands for broader bandwidth; the sources described emit GWs for weeks to years across the multiband to accumulate high signal-to-noise ratios. This suggests the possibility that joint detection would greatly improve the parameter estimation of the binaries. We examine B-DECIGO's ability to measure binary parameters and assess to what extent multiband analysis could improve such measurement. Using non-precessing post-Newtonian waveforms with the Fisher matrix approach, we find for systems like GW150914 and GW170817 that B-DECIGO can measure the mass ratio to within < 0.1%, the individual black-hole spins to within < 10%, and the coalescence time to within < 5 s about a week before alerting aLIGO and electromagnetic facilities. Prior information from B-DECIGO for aLIGO can further reduce the uncertainty in the measurement of, e.g., certain neutron star tidally-induced deformations by factor of ∼ 6, and potentially determine the spin-induced neutron star quadrupole moment. Joint LISA and B-DECIGO measurement will also be able to recover the masses and spins of intermediate-mass binary black holes at percent-level precision. However, there will be a large systematic bias in these results due to post-Newtonian approximation of exact GW signals. on their dipole radiation [46], the uncertainties in parameter estimation, and tests of general relativity [47]. Other proposals for multiband GW astronomy include those in Refs. [30,[48][49][50][51]. Nakamura et al. [22] have initiated examination of the precision with which the binary parameters can be determined by B-DECIGO, considering GW150914-like non-spinning BBHs. However, the measurability of BH spins, GW170817-like BNSs and the prospects for the multiband analysis with B- DECIGO have not yet been fully revealed. We shall assess how precisely we are able to measure the parameters of BNS and BBH inspirals with B-DECIGO, and explore how multiband B-DECIGO and aLIGO/ET (or LISA) measurement improves their parameter estimation and science cases over those using only single-band detection. Some previous studies of parameter estimation in the DECIGO mission can be found in Refs. [30,44,48,[52][53][54].
Assuming quasi-circular inspiraling binaries, the coalescing time t c , the instantaneous number of GW cycles N c ≡ f 2 (df /dt) −1 with the GW frequency f , and the dimensionless characteristic strain amplitude h c of the binary are estimated as (normalized to GW150914-like BBHs; see Sec. 3) [55,56] where D L is the luminosity distance between the observer and the binary, and the chirp mass M as well as GW frequency f are measured at the observer, accounting for cosmological effects. When binaries are at a cosmological distance, in the geometrical units, all mass scales are redshifted by a Doppler factor of 1 + z with the source's cosmological redshift z. As a result, for instance, the GW frequency f S and component masses m S i at the source location are related to those at the observer f O and m O i via f O = f S (1 + z) −1 and m O i = m S i (1 + z), respectively. In the rest of this paper, we always quote physical quantities measured at the observer, dropping the labels 'S' and 'O' (unless otherwise specified) 4 .
Equations (1) - (3) show that in the B-DECIGO band GW150914 and GW170817 were visible for ∼ 10 days and ∼ 7 yrs prior to coalescence with large numbers of GW cycles of  [56] completed by GW150914 and GW170817 (in the inspiral phase) are also depicted as references. Right: The detectable luminosity distance R L (m) for equal-mass inspirals as a function of their redshifted total masses m. We assume a four-year mission lifetime for B-DECIGO and LISA [2] and a detection SNR threshold at 8, for which we average over all-sky positions and binary orientation (see Eq. (12)); R L (m) with m fixed becomes smaller by a factor of √ 4ν for unequal-mass systems and larger by a factor of 2.5 for the optimal geometry. The luminosity distances of BBHs and a BNS observed by aLIGO and Virgo are also marked. These figures compliment the similar Figs. 2 and 3 in Ref. [22]. ∼ 10 5 and ∼ 10 7 , respectively. More importantly, both of their characteristic strains in this deci-Hz band are h c ∼ 10 −21 , which are well above the target dimensionless noise amplitude of B-DECIGO ∼ 10 −23 around 1 Hz [22]. This allows B-DECIGO to observe GW150914-and GW170817-like binary inspirals out to ∼ 60 Gpc (z ∼ 6) and ∼ 1 Gpc (z ∼ 0.2), respectively, assuming a detection (sky and polarization averaged) SNR threshold at 8 for a 4 yr mission lifetime 5 .
Intermediate-mass BBHs with the redshifted total mass m ∼ 10 4 M ⊙ can also stay in the B-DECIGO band for ∼ 1 d with ∼ 10 3 GW cycles before their final merger. Indeed, the GW frequency emitted at the innermost stable circular orbit (ISCO) of a Schwarzschild BH with (redshifted) mass M is

Parameter estimation with B-DECIGO: Multiband measurement
In matched filtering analysis, the measurement errors of binary parameters are classified into two categories; the statistical error due to the random noise in the detectors, and the systematic errors resulting from, e.g., inaccurate waveform modeling of the expected GW signals. To underline what could be expected from B-DECIGO measurements, next we consider the inspiral phase of aligned-spin BNSs and BBHs in quasi-circular orbits, which is approximated well by the post-Newtonian (PN) waveforms [59]. For B-DECIGO, the achievable precision of the statistical error can be very high, because the statistical errors scale as ∝ SNR −1 (1 + O(SNR −1 )) [60,61], and FIG. 1 shows that BNSs and BBHs will be visible by B-DECIGO with high SNRs 6 . At the same time, recall that the precision in statistical errors is determined by a combination of the SNR and the bandwidth over which a detector accumulates the SNR (see, e.g., Ref. [62]). Indeed, the high SNR of systems like GW150914 and GW170817 in the B-DECIGO band comes from the large numbers of GW cycles of ∼ 10 6 (recall Eq. (2)) accumulated in the much greater bandwidth during inspiral than that in the aLIGO band. This suggests the interesting possibility that we may be able to precisely measure, for instance, individual BH spins in GW150914-like BBHs in the B-DECIGO band, which can be hard to measure in the aLIGO band because of their strong degeneracy in the parameter dependence of the PN waveform [63][64][65][66][67][68]. The effects of BH spins come into the waveform at higher frequencies, but the broader bandwidth available in the lower B-DECIGO band might allow for tighter constraints.
The multiband analysis with B-DECIGO will further reduce the statistical uncertainties using only single-band analysis. The key point is that information from earlier B-DECIGO (or LISA) analysis can constrain the prior on the aLIGO (or B-DECIGO) analysis most naturally. This prior information refines the estimation of, for instance, the NS spin-and tidally induced deformations in BNS waveforms, from which we can infer the NS internal structure (equation of state) [69,70]. Once again, their precision is limited by a partial degeneracy between these effects and the mass ratios, as well as spins in the PN waveform [71][72][73][74][75][76]. If the mass ratio is already precisely constrained from the early inspiral in the lower B-DECIGO band, the measurement of such matter imprints in the waveform could be more precise in the higher aLIGO band.
Nonetheless, these improvement in statistical errors could be hampered by the systematic errors in the PN modeling of BNS and BBH inspirals, which is still an approximation of true general-relativistic (numerical-relativity) waveforms. Here it is important to recognize that the systematic measurement error is SNR independent [77] while the statistical measurement error scales with the inverse of SNR. Because the BNSs and BBHs observable by B-DECIGO will have high SNR, ultimately, the systematic mismodeling errors might be the limiting factor in their measurement.
The remainder of the paper quantifies the statistical and systematic parameter estimation errors in inspiraling BNSs and BBHs, using B-DECIGO and multiband measurements. We begin in Sec. 2 with a review of our methodology, including our PN waveform model and the Fisher matrix formalism for parameter estimation of GW signals. In Sec. 3, we present our main results for the statistical and systematic errors, respectively. We finish with two scientific cases in Sec. 4 that could be done in multiband GW astronomy with B-DECIGO: the redshift measurement of cosmological BNSs with GW observation alone, and the characterization of final remnant BHs using BBH inspirals.

Parameter estimation using post-Newtonian waveforms
As a first step toward the full problem of parameter estimation, for simplicity, we use as our GW signal model the up-to-date inspiral-only PN waveform with tidal (finite-size) corrections [59,71,78] neglecting merger and ringdown, and employ the semi-analytic Fisher information matrix formalism [63,77,79] for the signal analysis. We will also neglect the orbital motion of B-DECIGO and LISA, and average the signals over the all-sky positions and binary orientations [79,80]. Our analysis is limited to the aligned-spin inspirals in quasi-circular orbits because these are the only configurations for which the NS and BH tidal influences on the GW phase are computed in the PN approximation.

Tidally corrected non-precessing 3.5PN waveform
Motivated by the fact that matched filtering is more sensitive to the phase of the signal than its amplitude, we work with the frequency-domain "restricted" stationary phase approximation to the PN waveform, in which both higher-multipole components and PN corrections to the wave amplitude are ignored. After averaging over the all-sky positions and binary orientations, the resultant waveform reads [55] with the "Newtonian" amplitude scaled by the factor 2/5. The GW phase Ψ(f ) is the sum of two contributions: (i) spinning point-particle terms that are independent of the nature of the NSs or BHs comprising the binary, superposed with (ii) finite-size terms (depending on the nature of the NS or BH comprising the binary) that arise from the rotational deformation of an axially symmetric NS 7 , and the tidal response of the BH or NS in the binary on the other companion. For non-precessing BNSs, the GW phase Ψ(f ) may be expressed as where t c and Ψ c are the coalescence time and phase, respectively, and we introduce the orbital velocity v ≡ (πmf ) 1/3 : an O(v 2n ) term is of nth PN order. The first term, ∆Ψ pp 3.5PN , is the spin-independent, point-particle contribution up to 3.5PN order, derived in Ref. [62] 8 . The second term, ∆Ψ pp−spin 3.5PN , is the 3.5PN spin-dependent, point-particle contribution that includes linear spin-orbit effects [59,83], quadratic-in-spin effects [84], and cubic-in-spin 7 Following the tradition in PN waveforms, we have classified the corrections from BH rotational deformation here as "spinning point-particle" terms [81,82]. 8 In our analysis, we break with tradition and keep only terms ∝ ln(v) at 2.5PN order. The terms ∝ v 5 become constant in Ψ BNS (f ) due to cancellation with the overall factor of v −5 , and can be absorbed into Ψ c [62]. The same is applied with the 2.5PN terms in Ψ BBH (f ) below.

6/22
effects [85]. Using the dimensionless spin parameter χ i ≡ S S i · ℓ S /(m S i ) 2 defined in terms of the source-frame individual body's mass m S i and spin vectors S S i as well as the unit normal ℓ S to the orbital plane, ∆Ψ pp−spin 3.5PN is given as a function of v and (ν, χ s , χ a ) where χ s ≡ (χ 1 + χ 2 )/2 and χ a ≡ (χ 1 − χ 2 )/2. Note that positive (negative) values of χ i correspond to the aligned (anti-aligned) configurations with respect to the orbital angular momentum of the binary. Their explicit expressions are computed in Refs. [81,82].
The third term, ∆Ψ NS−QM 3.5PN , is the finite-size correction due to the rotational deformation of a NS. Restricted to the dominant effect, this is well characterized by the (dimensionless) NS quadrupole parameter The spin-induced quadrupole moment scalar Q [89] is fixed when the source-frame NS mass m S and equation of state are given, encoding the NS internal structure. Such spin-induced quadrupole-moment corrections to the GW phase start from 2PN order beyond the lowest PN term in Eq. (6), and we include them to the 3.5PN order [90] (assuming m 1 < m 2 ) 10 : where we define the "combined" dimensionless quadrupole parameters scaling as the square of the NS spins bỹ These parameters are conveniently chosen such that (a) the leading-order correction at 2PN order depends only onQ, (b)Q = 0 = δQ for a BBH because a spinning BH has κ i = 1 [84], and (c) It is also important to recognize that the parameters (Q, δQ) implicitly include the redshift factor (1 + z) 3 when using the NS masses at the observer m O i [75]. The last term, ∆Ψ NS−tidal

6PN
, is the finite-size correction due to the quadrupole tidal response of a NS. Restricted to slowly changing tidal fields, this response can be characterized by the (dimensionless) NS tidal deformability [71] (see also Ref. [92]). Similar to the quadrupole moment scalar Q, both the second (electric-type) Love number k 2 [93] and the source-frame NS radius R S are fixed when m S and the equation of state are given. Such tidal correction to the GW phase starts from 5PN order beyond the lowest PN term in Eq. (6), and we are only concerned with the leading-order (5PN) and next-toleading-order contributions (6PN) [94,95] 11 : where "combined" dimensionless tidal deformabilitiesΛ and δΛ are given by (once again assuming m 1 < m 2 ) [13, 72, 98] They have the convenient propertiesΛ(Λ 1 = Λ 2 = Λ) = Λ and δΛ(Λ 1 = Λ 2 = Λ) = 0 in the equal-mass limit ν = 1/4. The parametersΛ and δΛ in Eq. (9) are redefined to include the Doppler factor (1 + z) 5 for Λ 1,2 because we have used the NS masses at the observer m O i . We also remark that each NS tidal deformability can be related to the NS's spin-induced quadrupole-moment parameter using quasi-universal relations [99].
Meanwhile, the GW phase Ψ(f ) for aligned-spin BBHs may have the form The spin-(in)dependent point-particle terms ∆Ψ pp 3.5PN and ∆Ψ pp−spin are the same as those for a BNS in Eq. (6). The third term, ∆Ψ BH−tidal

3.5PN
, is the finite-size correction due to the tidal response of a BH. Restricted to slowly changing tidal fields, each BH in a BBH is tidally heated and torqued by its companion [78,100,101] 12 . These tidal contributions to the GW phase first appear at 2.5PN order for aligned-spin BBHs, and we keep the leading-order (2.5PN) and the next-to-leading-order (3.5PN) contributions [102], including the 2.5PN and 3.5PN contributions due to the energy and angular-momentum fluxes across the BH horizon, and the 3.5PN secular corrections to the binary's binding energy and GW luminosity (energy flux emitted to infinity) accumulated over the inspiral timescale. ∆Ψ BH−tidal 3.5PN is also the function of v and (ν, χ s , χ a ), which thus adds extra spin-dependent, finite-size contributions to Ψ BBH (f ). Its explicit expression is derived in Ref. [103]; note that the tidal heating and torquing for non-spinning BBHs start only from 4PN order, yielding ∆Ψ BH−tidal

Parameter estimation: Statistical errors
The parameter errors due to the overall effect of detector noise now have a firm statistical foundation (see, e.g., Refs. [63,79]). We assume that the GW signal observed in a detector (the so-called "template") is modeled by the sky-averaged 3.5PN waveformh(f ; θ) (see Eq. (5) with the set of physical parameters θ. We also assume that the noise in a detector is additive, stationary, Gaussian with zero means, and uncorrected with each other when considering a multiband network of GW detectors. 11 The spin-tidal coupling term starts at 6.5PN [96,97], which is negligible here. 12 The tidal Love numbers of slowly spinning BHs are all zero [93,104].

8/22
We begin with the single detector configuration. In the matched filtering analysis, the SNR (corresponding to the maximum correlation with the optimal filter) for the given timedomain signal h is defined by The bracket denotes the inner product weighted by the noise power spectrum density S h (f ) (asterisk " * " is used for complex conjugation) [79] ( where [f in , f end ] is the frequency range determined by the detector setup and property of signals; see Sec. 2.4. Note that the SNR of Eq. (13) automatically gives the averaged SNR over all-sky position and binary orientation [63,80], because of the "sky-averaging" factor of 2/5 i n the waveform of Eq. (5) 13 . Equation (12) can be recast in terms of D L to describe the observable range for a fixed ρ ave [62,105]. This was used to plot the right panel of FIG. 1.
For Gaussian noise and high-SNR sources (together with caveats [60,61]), the standard Fisher matrix formalism allows us to estimate the statistical errors δθ ≡ θ − θ 0 associated with the measurement, where θ and θ 0 are the best-fit parameters in the presence of some realization of noise and the "true value" of the physical parameters, respectively. In the high-SNR limit, δθ has a Gaussian probability distribution [60] where p (0) (θ) are the prior probabilities of the physical parameters; summation over repeated indices is understood (and we do not distinguish upper indices from lower ones). Here Γ ab ≡ (∂h/∂θ a |∂h/∂θ b )| θ=θ0 , is the Fisher information matrix defined in terms of Eq. (13), and its inverse defines the variance-covariance matrix Σ ab ≡ (Γ ab ) −1 for the Gaussian distribution of Eq. (14). Then, the root-mean-square error and cross correlations of parameters θ are given by where angle brackets denote an average over the Gaussian distribution of Eq. (14) (there is no summation over repeated indices here). By definition, each c ab must be restricted to the interval [−1, 1]; when |c ab | ∼ 1 (|c ab | ∼ 0), it is said that the two parameters θ a and θ b are strongly correlated (almost uncorrelated). Now we return to a multiband network configuration of detectors (e.g., "B-DECIGO + aLIGO", and so on). Because we assumed that the noise in the different detectors is uncorrelated, the total network SNR and Fisher information matrix are simply the sum of the 13 Because all "sky-averaging" factors are included in the waveform of Eq. (5), our convention for S h (f ) is the standard "non sky-averaging": where δ is a delta function,ñ(f ) is the Fourier component of the noise n(t), and angle brackets mean ensemble averaging with respect to the noise distribution.

9/22
individual (averaged) SNRs and Fisher information matrix for each detector: The total variance-covariance matrix for Eq. (14) is then given by Σ ab ≡ (Γ tot ab ) −1 , from which we can estimate the corresponding total root-mean-square error and cross correlations of parameters, making use of Eq. (15). Equations (14) and (16) directly show the advantage of parameter estimation with the multiband GW network. Having a priori knowledge from detector I in a different GW band, the parameter estimation with detector II could be more precise than a single-band analysis using only detector II.

Parameter estimation: Systematic errors
Next, we collect a few key results from the theory of GW signal analysis to measure the systematic mismodeling error; this arises from the fact that our PN waveform of Eq. (5) used in the statistical analysis only approximates the true general-relativistic signals.
We focus only on the waveform phasing error due to the neglect of the 4PN non-spinning point-particle term in the test-mass limit (ν = 0), which can over-dominate the error budget in measurement of NS tidal effects in the aLIGO band [72,98,106]. With this setup, we model the "true" GW signal by the sky-averaged PN waveformh T (f ) ≡ Af −7/6 e iΨT(f ) , making use of Eq. (5), and the true GW phase is Ψ The ν-independent coefficient |c 4 | ∼ 3200 is computed in Ref. [107] built on the results in Refs. [108,109]; the calculation of the ν-dependent correction to ∆Ψ pp 4PN is a current frontier in PN modeling [110][111][112].
A standard data-analysis-motivated figure of merit is the match [113,114] that measures the accuracy of the approximate (3.5PN) waveformh =h BNS/BBH of Eq. (5) by comparing to the true waveformh T with identical (true) source parameters θ 0 : where Here, maximizing over the phase ∆Ψ c is done analytically [115]. Waveform models with low match ( 0.97) are generally considered to be not "faithful" in the parameter estimation [116].
Still, how much systematic bias on the parameter estimation does the high-match ( 0.97) waveform generate? Given a best-fit waveformh(θ) to the detector output with the best-fit parameters θ, Cutler and Vallisneri showed that the systematic estimation errors in the source parameter ∆θ ≡ θ − θ 0 can be estimated by minimizing the inner product (ĥ T (θ 0 ) −ĥ(θ)|ĥ T (θ 0 ) −ĥ(θ)). In the high-SNR regime, it was shown that the minimization of this inner product yields [77] (see also Ref. [72]) ∆θ a = 3 32 where we use the sky-averaged PN waveform of Eq. (5) with the assumption ∆θ a ∂ a Ψ 1. In contrast to the statistical errors σ ∼ SNR −1 , it is important to recognize that the systematic 10/22 error ∆θ a is essentially independent of SNR because both the SNR of Eq. (12) and the Fisher information matrix scale as SNR 2 ∼ A 2 ∼ Γ ab , while ∆θ a ∼ A 0 .

Noise sensitivity of B-DECIGO and other GW detectors
The computation of the statistical and systematic parameter estimation errors will require as an input the model of noise power spectrum density S h (f ) corresponding to each GW detector. For B-DECIGO, we use the expected S h (f ) proposed by Nakamura et al. [22] (see also Ref. [53] for the (original) DECIGO configuration) with For other GW bands, we consider aLIGO [117] and a 3rd generation GW detector, the Einstein Telescope (ET) [118], in the hecto-Hz band as well as LISA [119] in the milli-Hz band (see also Ref. [120]). Their analytic fits to S h (f ) can be found in the cited references,

Binary parameters
In our simplified version of binary problems, the sky-averaged PN waveformsh BNS andh BBH in Eq. (5) depend on 11-and 7-dimensional parameters, respectively: θ BNS = (ln A, f 0 t c , Ψ c , ln m, ν, χ s , χ a ,Q, δQ,Λ, δΛ) , Here, we absorb all amplitude information into the single parameter A, and set f 0 = 1.65 Hz, at which S h (f ) of B-DECIGO is minimum. However, it is known that the PN waveform of Eq. (5) yields the block diagonal form of the Fisher information matrix Γ ln A a = δ ln A a ρ 2 ave , which means that ln A is entirely uncorrelated with the other parameters [64]. For this 14 The factor 20/3 = 5/ sin 2 (π/3) in Eq. (2.1) of Ref. [119] is a standard conversion factor between aLIGO and LISA configurations, which is (the inverse of) the products of 1/5 from an average over the pattern functions and 3/4 = sin 2 (π/3) from the angle between detector arms being 60 • [121]. Recall that our convention for S h (f ) is "non-sky-averaging" because the amplitude parameter A in Eq. (5) has already accounted for the pattern-average factor 1/5. To avoid counting this factor twice, we import S h (f ) from Ref. [119] after replacing the factor 20/3 with 4/3. 15 The abrupt cutoff of the waveform at f ISCO could artificially improve the parameter estimation if the waveform has sufficient noise-weighted power at the ISCO frequency [123]. For simplicity, we ignore this systematic bias.

11/22
reason, we shall remove ln A from our list of independent parameters and only consider the other 6-dimensional parameters, assuming that they are unconstrained (with "flat" priors p (0) (θ) ∼ const. in Eq. (14)) 16 .
In this work er consider the following three binaries for B-DECIGO; the parameter of System C is motivated by LISA's "threshold" system (with a source-frame total mass of m S ∼ 3000M ⊙ and mass ratio of m 1 /m 2 = 0.2) that defines LISA's observational requirement [2].
The coalescence time and phase are t c = 0.0 = Ψ c for each system. When we measure them using multiband network detectors, systems A, B, and C are also visible in the aLIGO band, aLIGO and LISA band, and LISA band, respectively (recall FIG. 1).

Results of parameter estimation
In what follows we summarize our parameter estimation results for each system. We stress that the methodology of our analysis is extremely simplified, e.g., with a flat prior and by using only the inspiral-only PN waveform; indeed, our aLIGO estimation errors for System A ("GW170817") and System B ("GW150914") obviously contradict those measured by aLIGO and advanced Virgo [8,13]. Thus, the estimation errors that we quote below should be only indicative and tentative. Our results have to be followed up by using more accurate inspiral-merger-ringdown waveforms with additional known effects (e.g., NS spins, spin-induced precession, etc.) and the orbital motion of B-DECIGO, a more rigorous Baysian-posterior-based parameter estimation method, and extended to a more exhaustive study of parameter spaces in future. 16 It should be borne in mind that in general a different prior assumption leads to different conclusions on the parameter estimation (see, e.g., Refs. [64,124]). From this point of view, it would be more physical to take into account priors for the fact that χ 1,2 and ν are restricted to the interval [−1, 1] and (0, 1/4], respectively. Analysis that assumes a certain prior forΛ and δΛ would also improve their estimation [13,73]. 17 For the GW170817-like equal-spin BNS (System A) described below, we have δQ/Q ∼ 0.07 and δΛ/Λ ∼ 0.08, and (δQ, δΛ) in the GW phase are further suppressed by the factor √ 1 − 4ν ∼ 0.04; they are essentially negligible in our analysis.

12/22
However, we are confident that the overall trend of our results (e.g., the order of magnitude of errors) should be robust, and provides a realistic idea of what can be measured with B-DECIGO as well as the multiband network including it.

Statistical errors: Fixed SNR
The statistical parameter estimation errors scale as σ ∝ SNR −1 (1 + O(SNR −1 )) [60,61], and their achievable precision depends on both the SNR and the bandwidth, over which the SNR is accumulated. For the former aspect, B-DECIGO measurement has already shown an advantage in FIG. 1, where we find that System A, B, and C are all high-SNR (∼ 10 2 ) sources. Furthermore, the multiband measurement does better than B-DECIGO alone because the systems considered are always much louder in the network SNR (16).
Here, we look at the latter aspects of the statistical errors, i.e., their improvement arising from the broader bandwidth. To best quantify this, we introduce the normalized statistical errors δθ normalized to a fixed SNR ρ, and display the statistical errors in terms of δθ; the overall statistical errors are simply recovered from σ = δθ/ρ. We start by looking at System A (GW170817-like BNS). In the left panel of FIG. 2, we show the B-DECIGO statistical errors δθ BD normalized to its SNR ρ BD ave (see Table 1) as a function of the lower cutoff frequency f in . Here we leave the errors δχ s,a and (δQ, δΛ) out of the plot because B-DECIGO is not able to discern the parameters associated with the waveforms; they are significant only in the higher aLIGO band, as anticipated. We see that δθ BD becomes significantly smaller as the signal contributes to broader bandwidth 18 . For a 4 yr observation before the coalescence, e.g., δν/ν is below 0.05% mark, and the overall statistical error in t c is σ tc 5.0 s when f end = 1.0 Hz (corresponding to ∼ 6 d before reaching at f ISCO ). Hence, B-DECIGO is able to alert aLIGO and electromagnetic observatories about the time of merger well in advance. Table 1 reports δθ for System A measured by two different multiband GW networks, "B-DECIGO + aLIGO" (the first column) and "B-DECIGO + ET" (the second column); they are now normalized to the corresponding total network SNR ρ tot . We see that δm and δν in the B-DECIGO measurements can be far better than those in aLIGO (ET) measurements by ∼ 2 orders of magnitude. For the multiband measurement with B-DECIGO + aLIGO (ET), the improvement in δm and δν is only incremental as B-DECIGO already measures them quite precisely during the early inspiral. When it comes to the NS (symmetric) spins χ s and matter imprints like quadrupole and tidal parameters (Q,Λ), however, the benefit from having B-DECIGO information becomes drastic. The multiband analysis 18 We find that δθ a BD in FIG. 2 can produce "bumps" when the corresponding covariance matrix c ab of Eq. (15) changes their signs, crossing zero. When c ab = 0, θ a and θ b become entirely uncorrelated with each other, and this suddenly "improves" (or "worsens") the parameter estimation. In our investigation, this should be another systematic bias due to PN waveforms because the varying PN order of the GW phase Ψ BNS/BBH in Eqs. (6) and (11) shifts the location of these bumps (for fixed binary parameters). We postpone a more detailed study of this issue to the future, but their systematic nature should lend caution to future parameter-estimation studies of binary inspirals based on the Fisher information matrix and PN waveforms. can reduce the error in the tidal parameter δΛ/Λ by an order of magnitude in aLIGO or ET alone. Furthermore, it is interesting to observe that there is the potential to extract the spin and spin-induced quadrupole parameters (χ s ,Q) from the GW waveforms making use of B-DECIGO + aLIGO (ET) measurements. While the errors (δχ s /χ s , δQ/Q) determined by only B-DECIGO or aLIGO(ET) are too large to yield any constraint on them, the multiband measurement significantly reduces these errors by two orders of magnitude. Indeed, these results benefit from the mass ratio being very precisely measured by B-DECIGO. Our estimation suggests that in the future B-DECIGO + aLIGO (ET) multiband era it will allow us to measure the NS equation of state much better relying on both spin and tidal effects.
We next turn to examine System B (GW150914-like BBH). In the right panel of FIG. 2 we plot δθ BD normalized to the B-DECIGO SNR ρ BD ave (see Table 2) as a function of f in . The overall trend of δθ BBH is similar as that of the left panel (System A), but, strikingly, B-DECIGO is able to measure the individual BH spins χ 1,2 ; after converting δχ BD s,a to the overall statistical errors of BH spins σ χ1,2 for ρ BD ave ) with the standard variance-propagation formula 19 , we find that χ 1,2 can be measured better than 10% (for a 4 yr observation). This result seems counterintuitive because the spin effects in the GW phase first enter through the 1.5 PN spin-orbit couplings [59], expected to emerge in the high frequency aLIGO band. However, recall that the spin-dependent point-particle contribution ∆Ψ pp 3.5PN in Eqs. (6) and (11) are proportional to the inverse of v −5 . Because of that, the 1.5PN spin-orbit term 19 Suppose random variables are collected in the random vector y, together with a non-linear function y = f (x) of Gaussian random variables x. In a first-order approximation of a Taylor series expansion of y = f (x 0 ) + Jdx . . . around a given vector x 0 with Jacobian matrix J ab = (∂f a (x))/∂x b | x=x0 , the variance-covariance matrix Σ yy for y is given by Σ yy = J Σ xx J T where Σ xx is the variance-covariance matrix for x.  Table 1 Statistical parameter errors δθ normalized to the total multiband SNR ρ tot for System A (GW170817-like BNS). The first and second columns give δθ for the two different multiband networks B-DCEIGO + aLIGO and B-DECIGO + ET, respectively (thus, using different total SNRs ρ tot for δθ). play an important role for B-DECIGO measurement from this point of view; its neglect adds ∼ 2% relative uncertainties when measuringχ s,a . Table 2 reports δθ for System B measured by three different multiband GW networks, "B-DECIGO + aLIGO" (the first column), "B-DECIGO + ET" (the second column), and "LISA + aLIGO" (the third column), each normalized to the total network SNR ρ tot . We see that the B-DECIGO measurement ofθ BBH can be all improved tremendously, by 3 to 4 (1 to 2) orders of magnitude, compared to the aLIGO (ET) measurement 20 . Althoughθ BBH is largely constrained in the B-DECIGO band, there is further several factors of improvement for joint B-DECIGO + aLIGO (ET) analysis. Also, the broader conclusions that we can draw from "LISA + aLIGO" measurement agree with those by Sasana [36] and Vitale [47], which focused on the similar GW150914-like system; LISA can predict the coalescence time within less than 30s, and prior information from LISA can reduce the uncertainties of individual BH spins in the aLIGO measurement (we shall not, however, attempt a precise comparison with them because of the significant difference in methodology). In principle, we could even consider full multiband measurement with a "LISA + B-DECIGO + aLIGO (ET)" network. However, it makes virtually no difference compared to the B-DECIGO + aLIGO (ET) network. B-DECIGO already measuresθ BBH very precisely, and the LISA contribution to ρ tot is very small. Finally, Table 3 reports δθ for System C (LISA's "threshold" BBH) measured by the multiband GW network "B-DECIGO + LISA", normalized to ρ tot . This clearly shows that B-DECIGO can measure θ BBH very well, except for χ a . We find that in terms of 20 Recall that they are overestimated due to our simplified inspiral-only analysis in the aLIGO band. If we instead compare our B-DECIGO results with the measurement uncertainties of GW150914 [125], the improvement against the aLIGO measurement may be by ∼ 2 orders of magnitude.   Table 3 Normalized statistical parameter errors δθ with respect to the total multiband SNR ρ tot for System C (LISA's "threshold" BBH). The notation is similar to that of Table 1. and they are evaluated at the true values, ν = 0.138889, χ s = 0.8, and χ a = 0.1.

Systematic errors
We note that for all of Systems A, B, and C, the inclusion of NS or BH spin effects to the GW phase is highly recommended although NS spins in waveform modeling are often neglected as they are much smaller than BH spins. Making use of the formula in Eq. (18), for example, the mismatch with and without the highest spinning point-particle terms at the 3.5PN order (Eq. (6b) of Ref. [82]) are all above the 3% mark when the body's spin is χ i 0.2; this value of spin is essentially independent of the choice of detectors, and becomes much smaller when not including the lower-order PN spin terms for the limit of 3% mismatch.
These results motivate the continued development of PN waveform modeling for accurate parameter extraction from future B-DECIGO and LISA measurement.

Discussion
For each system considered in this work, we have shown that multiband measurement with B-DECIGO will determine the binary parameters (masses, spins, NS quadrupole parameters and Love numbers) with percent-level precision, if the systematic bias due to PN waveform mismodeling is under control. Consequently, this exquisite precision will be able to enhance the already (or expected) rich payouts gained from aLIGO and LISA observation. Below we shall highlight two potential examples from cosmography and fundamental physics, focusing on the joint B-DECIGO + aLIGO (ET) measurement of the GW170817-like BNS and GW150914-like BBH, respectively.
The first example is the redshift measurement of GW170817-like BNS using only GW observation, proposed by Messenger and Read [132] as well as Harry and Hinderer [75]. Assuming that the NS equation of state is well constrained from the BNS waveform of the late inspiral and merger [76], the measurement of NS quadrupole parameters κ i and tidal deformabilities Λ i (i = 1, 2) imprinted in the (early) inspiral waveform allows determination of the source redshift z directly. Essentially, the point is that the parameters (κ i , Λ i ) manifestly depend on the source-frame masses m S i scaling as 5 , respectively, where m O i are the observer-frame NS masses. Because (κ i , Λ i ) are related with m S i through the NS equation of state, the measurement of (κ i , Λ i ) and m O i is then translated to m S i and m O i , from which we can determine the source redshift z.
This approach to cosmography benefits from the joint B-DECIGO + ET measurement. The Einstein Telescope will place strong constraints on the NS spins and tidal influence in BNSs, while B-DECIGO can precisely measure the NS masses: recall Table 1. Following the method of Ref. [132] to examine the uncertainties in our inspiral-only waveform parameters, we estimate θ BNS in Eq. (21) for System A relocated at ∼ 2.92 Gpc (z = 0.50), leaving (Q,Λ) out but with the addition of z, making use of Eqs. (15). The statistical errors for z then read δz/z = (1.17 × 10 −1 ) ET , (6.93 × 10 −2 ) ET+BD , where the SNRs are ρ ET ave = 6.62 and ρ tot = 7.00. Notably, we see that the uncertainty is ∼ 60% of what would be obtained with ET analysis alone. Joint B-DECIGO + ET measurement of BNSs would thus be quite valuable for this cosmography.
The second example is constraining the mass and spin of the final remnant BH after merger. For non-precessing BBHs, there are now numbers of numerical-relativity fitting formulas that enables the mapping of the initial BH masses m 1,2 and dimensionless spin parameters χ 1,2 to the final mass M f = M f (m, ν, χ s , χ a ) and dimensionless spin parameter χ f ≡ |S f |/M 2 f = χ f (m, ν, χ s , χ a ) of the remnant Kerr BH [133][134][135][136]. Because B-DECIGO (and ET) will measure m 1,2 and χ 1,2 very precisely, it will set stringent constraints on M f and χ f , making use of these fitting formulas.
(24) The joint B-DECIGO + ET analysis will further reduce the uncertainties in M f and χ f by up to factor of four in B-DECIGO + aLIGO analysis. Given that aLIGO and ET independently measure M f and χ f from the merger-ringdown phase, this drastic improvement in estimation from the inspiral phase will help to strengthen the inspiral-merger-ringdown consistency test of general relativity [139][140][141][142][143]. Other extreme gravity tests that can be done with binary inspirals, including direct and parameterized test of gravity theory [144][145][146] and the "nohair" test of BBH nature [90], might benefit from B-DECIGO and multiband measurement as well.
In conclusion, we expect B-DECIGO measurement of BNSs and BBHs in the deci-Hz band will complement, e.g., aLIGO, advanced Virgo, KAGRA, IndIGO, ET and LISA observation in the hecto-Hz and milli-Hz bands, boosting our understanding of astrophysics, cosmology, and gravity science. In the era of multiband GW astronomy, we shall decide and go, DECIGO [22,23].