Preliminary study on parameter estimation accuracy of supermassive black hole binary inspirals for TianQin

TL;DR

This study evaluates how accurately TianQin can estimate parameters of supermassive black hole binaries using Fisher matrix analysis, considering spin, eccentricity, and post-Newtonian effects, with simulations showing promising precision levels.

Contribution

It introduces a comprehensive Fisher matrix analysis including spin, eccentricity, and higher-order post-Newtonian effects for TianQin's parameter estimation of SMBH binaries.

Findings

Mass parameters can be estimated with sub-percent accuracy.

Luminosity distance can be measured within 1-3%.

Angular resolution is better than 12 square degrees.

Abstract

We use the Fisher information matrix method to calculate the parameter estimation accuracy of inspiraling supermassive black holes binaries for TianQin, a space-borne laser interferometric detector aimed at detecting gravitational waves in the millihertz frequency band. The `restricted' post-Newtonian waveform in which third order post-Newtonian (3PN) phase including spin effects (spin-orbit $\beta$ and spin-spin $\sigma$) and first-order eccentricity contribution is employed. Monte Carlo simulations using $10^3$ binaries for mass pairs with component masses in the range of $({10^5},{10^7}){M_ \odot }$ and cosmological redshift $z=0.5$ show that the medians of the root-mean-square error distributions for the chirp mass $M_c$ and symmetric mass ratio $\eta$ are in the range of $\sim 0.02\% - 0.7\% $ and $\sim 4\% - 8\% $, respectively. The luminosity distance $D_L$ can be determined to be $\sim 1\% - 3\% $, and the angular resolution of source $\Delta \Omega $ is better than 12 deg$^2$. The corresponding results for $z=1.0$ and $2.0$, which are deteriorated with the decreasing of the signal-to-noise ratio, have also been given. We show that adding spin parameters degrades measurement accuracy of the mass parameters (${M_c}$, $\eta$), and the time and the orbital phase of coalescence ($t_c$, $\phi _c$); the inclusion of the first-order eccentricity correction to the phase worsens the estimation accuracy comparing with the circular cases. We also show the effects of post-Newtonian order on parameter estimation accuracy by comparing the results based on second order and third order post-Newtonian phases. Moreover, we calculate the horizon distance of supermassive black hole binaries for TianQin.

Tables

Table 1. Table 1: The medians of the SNR distributions and the rms error distributions of parameters for 13 mass pairs, each of which has 10 3 superscript 10 3 10^{3} trials of SMBHBs. Sources locate at redshift z = 0.5 𝑧 0.5 z=0.5 , 1.0 1.0 1.0 , and 2.0 2.0 2.0 , respectively, with samples of inclination angle cos ⁡ ι ∼ 𝒰 ​ ( − 1 , 1 ) similar-to 𝜄 𝒰 1 1 \cos\iota\sim\mathcal{U}(-1,1) and angular positions of the sources cos ⁡ θ ∼ 𝒰 ​ ( − 1 , 1 ) similar-to 𝜃 𝒰 1 1 \cos{\theta}\sim\mathcal{U}(-1,1) and ϕ ∼ 𝒰 ​ ( 0 , 2 ​ π ) similar-to italic-ϕ 𝒰 0 2 𝜋 \phi\sim\mathcal{U}(0,2\pi) . The 3PN waveform that includes both spin and eccentricity effects is adopted. We set the eccentricity e 0 = 0.2 subscript 𝑒 0 0.2 {e_{0}}=0.2 at f 0 = 10 − 4 subscript 𝑓 0 superscript 10 4 {f_{0}}=10^{-4} Hz, t c = 0 subscript 𝑡 𝑐 0 t_{c}=0 sec, ϕ c = 0 subscript italic-ϕ 𝑐 0 \phi_{c}=0 rad, and β = σ = 0 𝛽 𝜎 0 \beta=\sigma=0 . The observation time T obs = 3 subscript 𝑇 obs 3 T_{\rm obs}=3 months.

($m_1,m_2$)	z	SNR	$\Delta \ln M_c$	$\Delta \ln D_L$	$\Delta \ln \eta$	$\Delta t_c$	$\Delta {\varphi }_c$	$\Delta \beta$	$\Delta \sigma$	$\Delta e_0$	$\Delta \Omega$
$\left(M_{\odot }\right)$			(%)	(%)	(%)	$\left(s\right)$	(rad)			$\left({10}^{-4}\right)$	$\left({10}^{-5}\mathrm{str}\right)$
$({10}^5,{10}^5)$	0.5	879	0.02	0.96	7.74	5.72	4.38	0.22	1.48	2.02	35
	1.0	437	0.04	2.25	15.4	14.6	8.73	0.44	2.96	1.68	205
	2.0	198	0.11	5.20	31.9	39.1	18.6	0.98	6.34	2.85	1062
$({10}^5,3\times {10}^5)$	0.5	917	0.03	0.99	5.51	8.22	4.56	0.16	1.19	1.49	36
	1.0	375	0.06	2.46	11.2	22.5	9.33	0.33	2.47	1.16	250
	2.0	135	0.19	5.50	22.2	73.5	19.7	0.84	5.58	1.71	1246
$(3\times {10}^5,3\times {10}^5)$	0.5	966	0.04	1.03	5.38	10.3	3.32	0.20	1.21	0.80	43
	1.0	369	0.08	2.50	10.5	27.1	6.53	0.42	2.43	0.65	238
	2.0	126	0.28	5.76	23.2	98.8	15.9	1.16	6.28	2.73	1304
$(3\times {10}^5,6\times {10}^5)$	0.5	781	0.06	1.09	4.66	15.7	3.56	0.22	1.23	0.50	46
	1.0	283	0.12	2.66	9.27	45.3	7.61	0.48	2.66	1.16	275
	2.0	102	0.41	6.86	20.8	172	19.0	1.38	7.00	5.12	1948
$(6\times {10}^5,6\times {10}^5)$	0.5	650	0.08	1.27	4.52	21.7	3.42	0.27	1.42	0.57	54
	1.0	274	0.16	2.68	9.26	58.3	7.11	0.57	2.89	1.91	309
	2.0	90	0.56	7.08	20.6	244	18.6	1.74	8.13	8.05	2064
$(6\times {10}^5,{10}^6)$	0.5	583	0.10	1.15	4.35	29.9	3.80	0.31	1.51	0.99	57
	1.0	197	0.22	2.87	8.65	93.5	8.15	0.70	3.29	2.93	325
	2.0	72	0.72	7.46	19.8	365	20.7	1.99	8.72	10.9	2374
$({10}^6,{10}^6)$	0.5	529	0.13	1.18	4.34	40.1	3.83	0.37	1.68	1.50	57
	1.0	203	0.26	2.96	8.64	113	7.97	0.76	3.48	3.82	361
	2.0	73	0.86	7.75	20.2	460	20.8	2.29	9.65	13.3	2377
$({10}^6,3\times {10}^6)$	0.5	310	0.23	1.63	3.82	123	5.74	0.54	2.08	3.27	107
	1.0	123	0.45	4.04	7.95	348	11.9	1.12	4.25	6.97	712
	2.0	39	1.27	11.3	16.9	1439	28.9	3.26	10.9	20.4	5485
$(3\times {10}^6,3\times {10}^6)$	0.5	274	0.29	1.69	4.13	185	4.98	0.70	2.49	4.53	120
	1.0	103	0.57	4.34	8.08	513	10.1	1.47	5.08	9.12	756
	2.0	35	1.86	12.4	19.6	2296	27.7	4.60	14.5	30.5	7708
$(3\times {10}^6,6\times {10}^6)$	0.5	224	0.38	2.08	3.82	336	5.90	0.87	2.78	6.08	205
	1.0	80	0.73	5.23	7.75	1058	12.6	1.90	5.94	12.1	1251
	2.0	27	2.54	13.9	20.3	4464	33.4	5.51	16.4	41.8	9742
$(6\times {10}^6,6\times {10}^6)$	0.5	191	0.46	2.22	3.94	493	5.76	1.02	3.13	7.45	246
	1.0	70	1.04	5.80	8.76	1454	13.2	2.46	7.23	17.1	1334
	2.0	23	3.46	16.2	21.8	7106	32.8	6.66	18.5	56.8	12144
$(6\times {10}^6,{10}^7)$	0.5	159	0.57	2.48	4.22	827	6.76	1.21	3.51	9.31	261
	1.0	55	1.45	7.04	9.88	2559	13.7	2.52	7.24	23.9	2274
	2.0	18	2.71	19.0	15.3	11016	38.6	7.88	20.8	45.1	18247
$({10}^7,{10}^7)$	0.5	137	0.73	2.75	4.81	1161	7.07	1.42	3.98	11.9	368
	1.0	50	1.46	7.88	8.33	3435	14.4	3.01	8.18	24.2	2830
	2.0	16	3.21	23.6	16.0	15135	44.4	10.5	26.8	53.5	28219

Table 2. Table 2: The medians of the rms error distributions (diagonal elements) and the mean absolute values of the correlation coefficients (off-diagonal elements) of parameters for Monte Carlo simulation of 10 3 superscript 10 3 10^{3} trials of 10 6 + 10 6 ​ M ⊙ superscript 10 6 superscript 10 6 subscript 𝑀 direct-product {10^{6}+10^{6}}{M_{\odot}} SMBHBs located at redshift z = 0.5 𝑧 0.5 z=0.5 . 3PN waveform is adopted. The rest of the setups are as in Table  1 . The table is symmetric; hence, we only list the upper triangular elements.

	$\Delta \ln M_c$	$\Delta \ln D_L$	$\Delta \ln \eta$	$\Delta t_c$	$\Delta {\varphi }_c$	$\Delta \theta$	$\Delta \varphi$
$\Delta \ln M_c$	0.011	0.229	0.907	0.522	0.684	0.242	0.273
$\Delta \ln D_L$	-	1.22	0.262	0.529	0.342	0.705	0.704
$\Delta \ln \eta$	-	-	0.143	0.595	0.807	0.273	0.310
$\Delta t_c$	-	-	-	8.06	0.663	0.488	0.612
$\Delta {\varphi }_c$	-	-	-	-	0.038	0.336	0.395
$\Delta \theta$	-	-	-	-	-	0.011	0.497
$\Delta \varphi$	-	-	-	-	-	-	0.015

Table 3. Table 3: As in Table  2 , except that the 3PN waveform including first-order eccentricity effect is adopted.

	$\Delta \ln M_c$	$\Delta \ln D_L$	$\Delta \ln \eta$	$\Delta t_c$	$\Delta {\varphi }_c$	$\Delta \theta$	$\Delta \varphi$	$\Delta e_0$
$\Delta \ln M_c$	0.028	0.255	0.950	0.623	0.761	0.265	0.316	0.425
$\Delta \ln D_L$	-	1.24	0.290	0.496	0.323	0.703	0.704	0.221
$\Delta \ln \eta$	-	-	0.243	0.686	0.851	0.303	0.357	0.545
$\Delta t_c$	-	-	-	10.1	0.719	0.455	0.593	0.435
$\Delta {\varphi }_c$	-	-	-	-	0.050	0.330	0.387	0.520
$\Delta \theta$	-	-	-	-	-	0.010	0.499	0.230
$\Delta \varphi$	-	-	-	-	-	-	0.016	0.271
$\Delta e_0$	-	-	-	-	-	-	-	0.195

Table 4. Table 4: As in Table  2 , except that the 3PN waveform including spin-orbit ( β 𝛽 \beta ) and spin-spin ( σ 𝜎 \sigma ) effects is adopted.

	$\Delta \ln M_c$	$\Delta \ln D_L$	$\Delta \ln \eta$	$\Delta t_c$	$\Delta {\varphi }_c$	$\Delta \theta$	$\Delta \varphi$	$\Delta \beta$	$\Delta \sigma$
$\Delta \ln M_c$	0.054	0.198	0.330	0.647	0.581	0.205	0.254	0.821	0.671
$\Delta \ln D_L$	-	1.26	0.192	0.343	0.208	0.702	0.692	0.203	0.212
$\Delta \ln \eta$	-	-	3.770	0.589	0.904	0.267	0.196	0.433	0.831
$\Delta t_c$	-	-	-	31.0	0.791	0.364	0.407	0.590	0.830
$\Delta {\varphi }_c$	-	-	-	-	2.812	0.270	0.242	0.543	0.986
$\Delta \theta$	-	-	-	-	-	0.012	0.506	0.228	0.264
$\Delta \varphi$	-	-	-	-	-	-	0.018	0.256	0.256
$\Delta \beta$	-	-	-	-	-	-	-	0.250	0.605
$\Delta \sigma$	-	-	-	-	-	-	-	-	1.131

Equations

$$h_{+}(t)$$

$$h_{\times }(t)$$

$$h(t)=F_{+}(\theta ,\varphi ,\psi )h_{+}(t)+F_{\times }(\theta ,\varphi ,\psi )h_{\times }(t).$$

$$F_{+}(t)$$

$$F_{\times }(t)$$

$$D_{+}(t;\theta ,\varphi )=$$

$$D_{\times }(t;\theta ,\varphi )=$$

$$D_L=\frac{1+z}{H_0}{\int }_0^z\frac{d{z}^{\prime }}{\sqrt{{\Omega }_M(1+z^{\prime }{)}^3+{\Omega }_{\Lambda }}},$$

$$M_c^{\prime }$$

$$M^{\prime }$$

$$\tilde{h}(f)=AQ{f}^{-7/6}{e}^{\mathrm{i}\Psi (f)},\mathrm{for}f>0,$$

$$A=-\sqrt{\frac{5}{96}}\frac{{M_c}^{5/6}}{{\pi }^{2/3}{D}_L},$$

$$Q=\sqrt{(1+\cos^{2}{\iota})^{2}F^{2}_{+}\big{(}t(f)\big{)}+(2\cos\iota)^{2}F^{2}_{\times}\big{(}t(f)\big{)}}\,.$$

$$t(f)=t_c-\frac{5}{256}{M_c}^{-5/3}{\left(\pi f\right)}^{-8/3}\sum _{k=0}^6{\tau }_k{x}^{k/2}+{\tau }_e,$$

$${\tau }_0=$$

$${\tau }_1=$$

$${\tau }_2=$$

$${\tau }_3=$$

$${\tau }_4=$$

$${\tau }_5=$$

$${\tau }_6=$$

$${\tau }_e=$$

$$\beta =$$

$$\sigma =$$

$$x=(\pi Mf{)}^{2/3}=(\pi M_{c}f{)}^{2/3}{\eta }^{-2/5}.$$

$$\Psi (f)=$$

$${\alpha }_0=$$

$${\alpha }_1=$$

$${\alpha }_2=$$

$${\alpha }_3=$$

$${\alpha }_4=$$

$${\alpha }_5=$$

$${\alpha }_6=$$

$${\varphi }_p(t(f))=\arctan \frac{-2\cos \iota F_{\times }(t(f))}{(1+{\cos }^2\iota )F_{+}(t(f))}.$$

$${\varphi }_D(t(f))=2\pi fR\sin \theta \cos \left(\frac{2\pi t(f)}{T}+{\varphi }_0-\varphi \right).$$

$${\varphi }_e(f)=-\frac{4239}{11696}(M_c\pi {)}^{-5/3}\frac{{f_0}^{19/9}}{f^{34/9}}{e_0}^2.$$

$$s(t)=h(t)+n(t)$$

$$(h_1|h_2)\equiv 4\mathrm{Re}{\int }_0^{\infty }\frac{{\tilde{h}}_1^{\ast }(f){\tilde{h}}_2(f)}{S_n(f)}df,$$

$$S_n(f)=\frac{S_x}{L_0^2}+\frac{4S_a}{(2\pi f{)}^4{L}_0^2}\left(1+\frac{10^{-4}\mathrm{Hz}}{f}\right),$$

Full text

I Introduction

The observational window of gravitational wave (GW) astronomy has been opened by the landmark detections by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo Abbott et al. (2016a, b, 2017a, 2017b, 2017c, 2017d); The LIGO Scientific Collaboration et al. (2018). Due to the limitations from gravity gradient noise and seismic noise, it is extremely challenging to detect GWs with frequencies $<10$ Hz by the terrestrial interferometers. In the low frequency band, the Laser Inteferometer Space Antenna (LISA) Amaro-Seoane et al. (2017) had been selected as the third large-class mission of the European Space Agency (ESA) with an anticipated launch around 2034. The key technologies for LISA, such as gravity reference system and space laser interferometry, have been successfully demonstrated by the LISA Pathfinder Armano et al. (2016).

TianQin is a proposed space-borne laser interferometric detector for gravitational waves in millihertz frequencies (0.1 mHz -1 Hz) Luo et al. (2016). The detector comprises three identical drag-free satellites orbiting around the Earth in a nearly equilateral triangular constellation. The geocentric distance of each satellite is approximately ${10^{5}}$ km which makes the arm length of the interferometer be about $1.73\times{10^{5}}$ km. The Keplerian orbit’s period of each satellite due to the gravitational attraction from the Earth is approximately $3.65$ days. Each pair of satellites is interconnected by two-way infrared laser beams which forms up to three Michelson interferometers. The guiding center of the constellation coincides with the geocenter and it moves around the Sun in the ecliptic orbit. The normal of the detector plane points toward the tentative reference source RX J0806+15, which is among the strongest GW sources in Galactic ultra-compact white-dwarf binaries Israel et al. (2002); Strohmayer (2008); Huang et al. (2019). The fundamentals of the satellite orbits and the response of each Michelson interferometer for TianQin have been studied in Hu et al. (2018). The nominal orbits and a set of alternative orbits have been optimized such that the stability requirements on arm length variation, relative velocity and breathing angle of the triangular constellation can be satisfied for a five-year mission lifetime Ye et al. .

As a millihertz frequency gravitational wave observatory, the design of the architecture for the TianQin mission and the trade-off among a variety of technologies for the instruments are driven by the attainable science objectives. The main categories of GW sources in TianQin’s frequency band Hu et al. (2017) are Galactic ultra-compact binaries, coalescing supermassive black hole binaries (SMBHBs), capture of stellar-mass compact objects by a massive black hole (MBH), i.e., extreme-mass-ratio inspirals (EMRIs) Hughes (2001), inspiral of stellar mass black hole binaries Sesana (2016), and stochastic GW background originating from primordial BHs Di and Gong (2018) and/or cosmic strings Ölmez et al. (2010). Among these, SMBHB mergers are arguably the most powerful GW sources; therefore, they deserve detailed investigations.

Observations indicate that the center of almost every galaxy hosts a SMBH whose mass is in the range of $10^{5}-10^{10}{M_{\odot}}$ Kormendy and Richstone (1995). In the hierarchical merger scenario of galaxy formation, large galaxies are assembled through multiple mergers during their lifetime. As they merge, the center SMBHs will approach each other due to the deep gravitational potential, dynamic friction, triple interaction, gravitational waves, etc. (Milosavljević and Merritt, 2003; Komossa, 2003). SMBHB can even be formed directly inside the first galaxies under certain physical conditions (Bromm and Loeb, 2003). SMBHB merger rates along the cosmic history predicted by three different population models and the science capability for different eLISA designs have been discussed Klein et al. (2016). The work on estimating SMBHB detection rates for TianQin based on the semi-analytical model is underway Wang et al. (2019). The science potential of TianQin on testing the black hole no-hair theorem by using the ringdown signals from SMBHB mergers has been studied Shi et al. (2019).

In addition to detection rates, it is important to forecast how accurately that TianQin can measure the parameters pertinent to SMBHB systems which potentially can be used to enable the subsequent studies, for example, formation and growth mechanism of seed black holes Madau and Rees (2001); Bromm and Loeb (2003), co-evolution of SMBHBs with host galaxies Magorrian et al. (1998); Kormendy and Ho (2013) and cosmography Hughes and Holz (2003); Tamanini et al. (2016). From the perspective of data analysis, one needs matched filtering to extract the deterministic signals, such as the ones from inspiral stage of SMBHBs, from the noisy data. This involves passing the detector’s strain data output through a bank of templates that are characterized by the signal parameter set (Kay, 1998; Helstrom, 1968). In general, Monte-Carlo simulation using synthetic data and a data-analysis pipeline will be required to systematically evaluate the performance of parameter estimation and signal reconstruction for a specific detector. However, this procedure is cumbersome and computationally demanding. For the ease of use and computational efficiency, the Fisher information matrix (FIM) method has been implemented to bound the parameter estimation accuracy of post-Newtonian (PN) inspirals for both stellar-mass black hole binaries detected by the ground-based detector LIGO Cutler and Flanagan (1994) and SMBHBs detectable by the space-borne detector LISA Cutler (1998). It has been shown that the results based on the FIM method are consistent with the ones from more sophisticated Bayesian parameter estimations for high signal-to-noise ratio (SNR) cases Rodriguez et al. (2013); Porter and Cornish (2015).

In the applications of FIM method, Cutler and Flanagan Cutler and Flanagan (1994) used the restricted 1.5PN inspiral waveform that includes spin-orbit parameter ($\beta$) to discuss the estimation accuracy of luminosity distance $D_{L}$, chirp mass $M_{c}$, reduced mass $\mu$ and $\beta$ for LIGO and Virgo network. By extending the waveform template phasing to 2PN and including both spin-orbit ($\beta$) and spin-spin ($\sigma$) parameters, Poisson and Will Poisson and Will (1995) found that the 1.5PN phasing actually underestimated the root-mean-square (rms) errors in $M_{c}$, $\mu$, and $\beta$. Furthermore, Arun et al. Arun et al. (2005); Arun (2006) adopted 3.5PN phasing, however, ignored spin and eccentricity effects, to emphasize the importance of employing higher order PN correction for parameter estimation.

For LISA, Cutler Cutler (1998) first used the waveform in Cutler and Flanagan (1994) to calculate the sky location error, which was later extended to 2PN waveform including spin effects by Berti et al. Berti et al. (2005). Lang et al. Lang and Hughes (2006); Lang et al. (2011) used 2PN waveform including the precession-induced modulations with partially aligned spins and higher PN harmonics to examine the impact upon parameter estimation for LISA. They found that the additional precession periodicity and higher PN harmonics can improve the accuracy. The first-order phase correction due to orbital eccentricity has been considered in Królak et al. (1995). Afterwards, it has been found that high eccentricity ($e_{0}\geq 0.6$) not only enhances the SNR but also improves LISA’s angular resolution of SMBHBs than the circular obits for relatively high equal mass systems ($\sim 10^{7}M_{\odot}+10^{7}M_{\odot}$) Mikóczi et al. (2012); Gondán et al. (2017).

In the current work, we mainly focus on the investigation of parameter estimation accuracy for SMBHB inspirals that can be achieved by TianQin. We employ the restricted post-Newtonian waveform that includes up to third order (3PN) phase and contributions from spin effects (both spin-orbit and spin-spin) and first-order eccentricity effects. The plan for the rest of the paper is as follows. In Sec. II, we give the GW waveform used in this work. In Sec. III, we briefly summarize the characteristics of the TianQin detector and its prospects for detecting GWs from inspiraling SMBHBs. In Sec. IV, we review some basics of the FIM method. Sec. V carries out the Monte-Carlo simulations for typical SMBHBs and presents the main results of the rms error distributions of estimated parameters and the comparison with variations in employed waveforms. Conclusions and discussions on possible future studies are given in Sec. VI. In Appendix A, we provide the details of the transforming GW signal in time-domain to frequency-domain based on the stationary phase approximation (SPA), which leads to validating the requirements of SPA for TianQin in Appendix B. Throughout this paper we adopt the geometrical units in which $G=c=1$.

II Waveform model

II.1 The time-domain waveform

In the source rest frame, we construct the center-of-mass coordinates $\{x_{1},x_{2},x_{3}\}$ and consider a binary consisting of two masses $m_{1}$ and $m_{2}$ in a circular orbit on the $x_{1}-x_{2}$ plane. We assume that there is no precession; thus, the orbital angular momentum vector ${\bf L}$ points along a fixed direction, i.e. the $x_{3}$ axis (here we adopt the conventions in Apostolatos et al. (1994); Creighton and Anderson (2011)). The source locates at a distance of $r$ with an orbital inclination angle $\iota$ defined as $\cos\iota=-\bf\hat{L}\cdot\bf\hat{N}$, where ${\bf\hat{N}}$ is the unit vector pointing towards the source from the detector and ${\bf\hat{L}}$ is the unit vector of ${\bf L}$. Using the Newtonian mass quadrupole formula, we can obtain the time-domain waveforms for the two polarizations of the GWs propagating along the ${\bf-\hat{N}}$ direction:

[TABLE]

where $M_{c}=\mu^{3/5}M^{2/5}=\eta^{3/5}M$ is the chirp mass (the total mass $M={m_{1}}+{m_{2}}$, the reduced mass $\mu={m_{1}}{m_{2}}/M$ and the symmetric mass ratio ${\eta=\mu/M}$), $t_{c}$ and ${\varphi_{c}}$ are the time and the orbital phase of coalescence, respectively.

For detecting continuous GW signals by a space-borne detector, it may be more convenient to work in the heliocentric ecliptic frame. The GW strain $h(t)$ depends on the detector’s antenna response to the two GW polarizations ${h_{+,\times}}(t)$,

[TABLE]

Here $F_{+,\times}$ are the antenna pattern functions. $(\theta,\phi)$ denotes the source’s ecliptic colatitude and longitude, $\psi$ is the polarization angle between one of the semi-major axes of the ellipse of GW polarization and the line of nodes. The orbit and orientation of the detector determine the antenna response to the incoming GWs. For the preliminary concept of TianQin Luo et al. (2016), the response of a Michelson-type interferometer that is valid in the full range of the interested frequencies has been given in Hu et al. (2018). Specifically, the antenna pattern functions can be written as follows:

[TABLE]

In the low-frequency region ($f<f_{\ast}\approx 0.28$ Hz) that is most relevant to the GWs from an inspiraling SMBHB, $F_{+,\times}$ becomes independent of the GW frequency $f$, and $D_{+,\times}$ can be analytically approximated as Hu et al. (2018)

[TABLE]

Here $\kappa_{1}=2\pi f_{\rm{sc}}t+\kappa_{0}$, $f_{\rm{sc}}=1/(3.65$ day), and $\kappa_{0}$ is the constant phase determined by the setup of the satellites’ coordinates (see Hu et al. (2018) for details). $(\bar{\theta}=1.65,\bar{\phi}=2.10)$ are the colatitude and longitude of the reference source RX J0806+15 in the heliocentric-ecliptic frame towards which the normal of the detector’s plane is pointed Luo et al. (2016).

As we will see in Sec. III, TianQin can detect SMBHBs located at cosmological distances; therefore, it is natural to replace the distance variable $r$ in Eq. (1) by luminosity distance $D_{L}$. In the standard flat $\Lambda$CDM cosmological model, $D_{L}$ can be expressed as a function of the cosmological redshift $z$ of the source as

[TABLE]

where the matter density $\Omega_{M}=0.32$, the dark energy density $\Omega_{\Lambda}=0.68$, and Hubble constant ${H_{0}}=67~{}{\rm km~{}s}^{-1}{\rm Mpc}^{-1}$ Calabrese et al. (2017). Correspondingly, the chirp mass and total mass measured in the source rest frame can be related to their redshifted counterparts by

[TABLE]

For simplicity, hereafter we will omit the prime symbol and redefine ${M_{c}}$ and $M$ as the redshifted chirp mass and total mass that are measured in the detector’s frame unless otherwise specified.

II.2 The frequency-domain waveform

In this subsection we compute the Fourier transform of the time-domain GW signal. We extend Eq. (1) to the “restricted” PN waveform for which the amplitude is kept to the leading order quadrupole term while the phase is kept to higher PN orders, since it is more important to make the phase coherent in GW signal detection.

Given a GW strain signal ${h(t)=A(t){\cos\Phi(t)}}$, we can obtain its Fourier transform $\tilde{h}(f)$ analytically using the stationary phase approximation (SPA, see Appendix A for details) under required constraints (see Appendix B for the validation) as follows,

[TABLE]

with $\mathrm{i}^{2}=-1$ and

[TABLE]

and

[TABLE]

Here and hereafter, we express $t$ as a function of $f$ in the frequency-domain waveform. For the phase of waveform that includes 3PN, spin and eccentricity effects, $t(f)$ can be explicitly written as Buonanno et al. (2009); Królak et al. (1995)

[TABLE]

with the coefficients

[TABLE]

Here, $\gamma_{E}=0.577$ is the Euler constant, $e_{0}$ is the eccentricity of the binary Keplerian orbit at the fiducial frequency ${f_{0}}$. Spin-orbit ($\beta$) and spin-spin ($\sigma$) parameters can be expressed as Berti et al. (2005); Blanchet et al. (1995)

[TABLE]

where $\boldsymbol{\chi}_{i}={\bf S}_{i}/m_{i}^{2}$ $(i=1,2)$ is the dimensionless spin parameter and ${\bf S}_{i}$ is the spin angular momentum for the $i$-th black hole. ${\bf\hat{L}}$ is the unit vector of the orbital angular momentum of the binary. For black holes $|\boldsymbol{\chi}|\leq 1$, $\left|\beta\right|\lesssim 9.4$ and $\left|\sigma\right|\lesssim 2.5$. As an example, if ${\bf S}_{i}$ aligns with ${\bf\hat{L}}$ and ${\chi_{i}}=0.2$, then $\beta=1.6$ and $\sigma=0.1$ for equal-mass binary. $x$ is the PN parameter

[TABLE]

The 3PN GW strain phase evolution that includes the polarization modulation, Doppler modulation and the first-order eccentricity correction is given by Buonanno et al. (2009); Cutler (1998); Królak et al. (1995)

[TABLE]

with the coefficients

[TABLE]

where $x_{0}$ is the PN parameter (Eq. 13) evaluated at the last stable orbit.

The polarization modulation is defined as

[TABLE]

The motion of the TianQin detector around the heliocenter will cause Doppler modulation on the GW phase. Following (Hu et al., 2018), it can be expressed as

[TABLE]

Here $\phi_{0}$ is the ecliptic longitude of the satellite guiding center at $t=0$, $T$ is one sidereal year and $R=1$ AU.

Although binary orbit can be effectively circularized by the time of coalescence in high frequencies for stellar-mass BHBs, the orbital eccentricity may not be negligible for SMBHBs during inspiraling in low frequencies. Potentially, the eccentricity can be detected by space-borne detectors. Here we adopt the first-order correction of GW phase due to eccentricity Królak et al. (1995):

[TABLE]

III Signal detection

The strain output of a detector,

[TABLE]

consists of a time series of the GW stain signal $h(t)$ and the detector’s equivalent strain noise $n(t)$. The noise is assumed to be sampled from a Gaussian stochastic process. We define the noise-weighted inner product $({h_{1}}|{h_{2}})$ of two data ${h_{1}}(t)$ and ${h_{2}}(t)$ as Finn (1992)

[TABLE]

where ${\tilde{h}_{1}}(f)$ and ${\tilde{h}_{2}}(f)$ are the Fourier transforms of ${h_{1}}(t)$ and ${h_{2}}(t)$. ${{S_{n}}(f)}$ is the one-sided power spectral density (PSD) of $n(t)$ for which the proposed functional form for TianQin is provided in Luo et al. (2016); Hu et al. (2018)

[TABLE]

where $L_{0}=1.73\times 10^{5}~{}{\rm{km}}$ is the arm length. $S_{x}=10^{-24}~{}{\rm{m}}^{2}/{\rm{Hz}}$ and $S_{a}=10^{-30}~{}{\rm{m}}^{2}{\rm{s}}^{-4}/{\rm{Hz}}$ are the PSDs of the position noise and residual acceleration noise, respectively.

The detectability of a given signal $h(t)$ is determined by the optimal signal-to-noise ratio (SNR) Finn (1992); Moore et al. (2015)

[TABLE]

where $h_{c}(f)=2f|{\tilde{h}}(f)|$ is the characteristic strain of the signal and $h_{n}(f)=\sqrt{fS_{n}(f)}$ is the counterpart of the noise. Note that both of these are dimensionless.

The average squared strain response for TianQin in Eq. (9) can be obtained by integrating over all possible sky locations ($\theta\in[0,\pi],\phi\in[0,2\pi]$), polarizations ($\psi\in[0,\pi]$), source orientations ($\iota\in[0,\pi]$), and the time within one orbital period of the satellites around the geocenter (${T_{\rm{sc}}}=3.15\times{10^{5}}~{}{\rm{s}}$):

[TABLE]

The frequency-domain waveform (see Eq. (7)) that adopts the average strain response $\langle|Q|^{2}\rangle^{1/2}$ is given by

[TABLE]

Combining Eqs. (21), (22) and (24), the SNR can be reformulated as

[TABLE]

where ${f_{\rm fin}}=\min({f_{\rm ISCO}},{f_{\rm end}})$ with the GW frequency at the innermost stable circular orbit ${f_{\rm ISCO}}={1/(6^{3/2}M\pi)}$ Hz and the upper cutoff frequency for TianQin ${f_{\rm end}}=1~{}{\rm{Hz}}$, ${f_{\rm in}}=\max({f_{\rm low}},{f_{\rm obs}})$ with the lower cutoff frequency ${f_{\rm low}}={10^{-5}}~{}{\rm{Hz}}$ and the initial observation frequency ${f_{\rm obs}}=4.15\times{10^{-5}}(M_{c}/10^{6}~{}M_{\odot})^{-5/8}(T_{\rm obs}/1~{}{\rm yr})^{-3/8}$ Hz. We will choose $T_{\rm obs}=3$ month which is the time window of each separated observation session. This is mainly due to the unique feature of TianQin’s mission operation which is intended to reduce the interference from the sunlight on the optical links and simplify the thermal control of the satellites (Luo et al., 2016).

Fig. 1 shows the characteristic strain of the detector noise $h_{n}(f)$ for TianQin (red line) as well as the most recent one for LISA (green line) Babak et al. (2017), and the characteristic strains of the signal $h_{c}(f)$ for three typical equal-mass (${10^{5}+10^{5}}{M_{\odot}}$, ${10^{6}+10^{6}}{M_{\odot}}$, ${10^{7}+10^{7}}{M_{\odot}}$) SMBHBs located at $z=0.5$. Here and hereafter, we refer the mass of SMBHB to the value measured in the source rest frame considering its direct connection with astrophysical investigations on the mass dependence of detection rates. However, in our calculation of SNR and FIM, we take into account the cosmological effects discussed in Sec. II.1. The blue lines with $h_{c}\propto f^{-1/6}$ present the GW contributions from the inspiral stage of the binaries which terminates at ${f_{\rm ISCO}}$. This is the most relevant part to the current work. Additionally, for the purpose of illustration and comparison, the corresponding ones from the merger (green lines) and ringdown (purple lines) stages have also been plotted based on the phenomenological waveform Ajith et al. (2011).

For $h_{c}\propto M_{c}^{5/6}\propto M^{5/6}$, the overall amplitude of the signal increases almost linearly with the total mass $M$, while the line that denotes the inspiral stage owns a lower cut-off frequency for a more massive system (${f_{\rm ISCO}}\propto M^{-1}$). Consequently, the competition between the integrand and the integration limits in Eq. (25) will result in a maximum of $\rho$ for a given $D_{L}$. For equal-mass systems, this is demonstrated in the upper panel of Fig. 2, in which four $D_{L}$ corresponding to $z=0.1$, $0.5$, $1.0$ and $2.0$ (from top to bottom) are used. The maxima of $\rho$ locate at $\sim$ a few $\times{10^{5}}{M_{\odot}}$.

The lower panel of Fig. 2 shows the horizon distance in terms of cosmological redshift (left axis) and luminosity distance (right axis) as a function of $M$ for fixed $\rho=10$, $50$, and 100 (from top to bottom), which may be regarded respectively, as criteria of the weak signal that can be marginally detected, the intermediate signal, and the strong signal. Specifically, for marginal detections with $\rho=10$, we can see that TianQin is capable of detecting sources as far as $z>30$ for the lower end of the mass spectrum of the SMBHB systems that are assembled in the early Universe. The detection space shown in Fig. 2 is largely overlapped with eLISA Hu et al. (2017); eLISA Consortium (2013).

Fig. 3 shows the fraction of the total SNR that can be accumulated in the observation time intervals before the merger of the equal-mass binary located at $z=1.0$. The total SNR is estimated by assuming TianQin’s continuous operation of 5 years (without the separation of $T_{\rm obs}=3$ months observation sessions). The minimum of the $99\%$ of the total SNR curve is determined by the interplay between the time to coalescence of a system $\tau$ ($\propto M^{-5/3}f^{-8/3}\eta^{-1}$) Vecchio (2004) and the characteristic strains of the detector noise $h_{n}(f)$ in Fig. 1. For $10^{6}+10^{6}{M_{\odot}}$ SMBHBs, one can see that more than $99\%$ of the total SNR can be obtained within one hour before the merger. Notice that over the vast majority of the mass range, the signal from the last few days’ observation will dominate the total SNR. Therefore, it is likely that, in practice, the signal from a SMBHB will either be detected by TianQin within $T_{\rm obs}$ or completely missed.

IV Fisher information matrix

Let ${\bar{\lambda}^{i}}$ be the true value of parameter ${\lambda^{i}}$ and $\hat{\lambda}^{i}={\bar{\lambda}^{i}}+\Delta{\lambda^{i}}$ be the estimated value from the data in the presence of noise. $\hat{\lambda}^{i}$ can be obtained from an estimator, such as the maximum-likelihood estimator (MLE) Kay (1993). Suppose that the signal is sufficiently strong; the joint probability density function for the estimation error $\Delta\boldsymbol{\lambda}$ can be approximated by a multivariate Gaussian distribution Cutler and Flanagan (1994); Creighton and Anderson (2011)

[TABLE]

where

[TABLE]

is the $(i,j)$ entry of the Fisher information matrix (FIM). Particularly, the rms error of parameter ${\lambda^{i}}$

[TABLE]

where the variance-covariance matrix (or simply the covariance matrix) $\Sigma={\Gamma^{-1}}$. The correlation coefficient between $\Delta\lambda^{i}$ and $\Delta\lambda^{j}$ can be defined as

[TABLE]

It is a dimensionless ratio that indicates the degree to which $\Delta\lambda^{i}$ and $\Delta\lambda^{j}$ are linearly correlated.

In this work, we focus mainly on the following parameters of the inspiraling SMBHB signal included in Eq. (7)

[TABLE]

The derivative parameters not presenting explicitly in Eq. (7) can be obtained by the propagation of errors. For example, the rms error associated with the reduced mass is

[TABLE]

In addition, the error in the sky localization $\Delta\Omega$ can be expressed in terms of the rms errors in $\theta$ and $\phi$ Berti et al. (2005)

[TABLE]

Note that both Eq. (31) and Eq. (32) take into account the correlations between parameters.

The derivatives of ${\tilde{h}}$ with respect to the parameters used in Eq. (27) are listed below:

[TABLE]

The explicit expressions for ${\partial\Psi}/{\partial{\lambda^{i}}}$ and ${\partial Q}/{\partial\lambda^{i}}$ on the right-hand sides of the above equations are too lengthy to show here.

In the following calculations, we neglect $\partial A/\partial\ln M_{c}$, $\partial Q/\partial\ln M_{c}$, and $\partial Q/\partial\ln\eta$, since these terms are ignorable comparing with the contributions from the terms related to $\Psi$. We will discuss this point with further details in Sec. VI.

V Simulations and results

In this section we use the FIM method to study the parameter estimation accuracy of inspiraling SMBHBs for TianQin. We assume that the strain data of the detector are from a Michelson-type configuration which consists of two-way optical links along two arms of the satellite constellation Hu et al. (2018). The functional form of instrumental noise PSD given in Eq. (21) and the associated parameter values are used for TianQin.

We carry out Monte Carlo simulations of $10^{3}$ SMBHBs for mass pairs listed in the first column of Table 1. The component masses are chosen such that the average SNR of those systems can be approximately larger than 20 (cf. Fig. 2). In this work, SMBHBs are located at $z=0.5$, $1.0$, and $2.0$ with the inclination angle $\cos\iota\sim\mathcal{U}(-1,1)$ and the angular positions of the sources $\cos{\theta}\sim\mathcal{U}(-1,1)$ and $\phi\sim\mathcal{U}(0,2\pi)$. Here, $\mathcal{U}(a,b)$ denotes the continuous uniform distribution between interval $[a,b]$. We set the eccentricity of the binary ${e_{0}}=0.2$ at ${f_{0}}=10^{-4}$ Hz. As in Cutler and Flanagan (1994), we are mainly concerned with the effects of the spin on the other parameters and choose $\beta=\sigma=0$ in all cases.

V.1 Parameter estimation of 3PN phase including spin and eccentricity effects

Table 1 shows the medians of the rms error distributions of the parameters employed in the 3PN waveform including both spin and eccentricity. Generally, the rms errors increase with the decreasing of the SNR for a specific detector. The peaks of the SNR can be read from the top panel in Fig. 2. One can see that the more distant sources would peak at less massive systems which merge at higher frequencies in the source frame and are subsequently redshifted to lower frequencies in the observer’s frame.

Specifically, for $z=0.5$ (corresponding to ${D_{L}}\approx 3~{}{\rm{Gpc}}$), we can see that $\Delta\ln M_{c}$ and $\Delta\ln\eta$ can be measured to be $\sim 0.02\%-0.7\%$ and $\sim 4\%-8\%$, respectively. It turns out that $\Delta\ln M_{c}$ is more sensitive to SNR than $\Delta\ln\eta$. $\Delta\ln D_{L}$ can be determined to be $\sim 1\%-3\%$. Despite that $D_{L}$ appears in the signal model as an overall factor, its estimation error does not simply scale as ${1/\rho}$. This is mainly due to the correlations between $\Delta\ln{D_{L}}$ and $\Delta\lambda_{i}$ (${\Sigma_{\ln{D_{L}}~{}\lambda_{i}}}\neq 0$ for $\lambda_{i}\neq\ln{D_{L}}$). Note that it is different from the overall factor ${\cal A}$ adopted in Berti et al. (2005); Cutler and Flanagan (1994); Poisson and Will (1995); Arun et al. (2005) which is entirely uncorrelated with the other parameters.

For $z=0.5$, the sky-position resolution $\Delta\Omega$, which is crucial for multi-messenger observations, can be measured to be $\sim 1-12~{}{\deg^{2}}$ ($10^{-5}~{}{\rm{str}}\approx 1/30~{}{\deg^{2}}$). Furthermore, the error for the estimated time of coalescence $\Delta t_{c}$ is less than 20 min for all cases, which may enable TianQin to send out prompt alerts to ground and space telescopes to search for the potential electromagnetic counterpart within $\Delta\Omega$.

For comparison, we calculated the SNR and parameter estimation accuracy for LISA based on its most recent noise curve Babak et al. (2017) (see Fig. 1). The medians of the SNR read 2267, 1358, and 180 for the three equal-mass (${10^{5}+10^{5}}{M_{\odot}}$, ${10^{6}+10^{6}}{M_{\odot}}$, ${10^{7}+10^{7}}{M_{\odot}}$) SMBHB systems located at $z=0.5$. As an example, the medians of the rms errors of $\Delta\ln M_{c}$, $\Delta\ln\eta$, $\Delta\ln D_{L}$, and $\Delta\Omega$ for the ${10^{5}+10^{5}}{M_{\odot}}$ system are $\sim 0.046\%$, $\sim 33\%$, $\sim 3.7\%$, and $\sim 26~{}{\deg^{2}}$, respectively. Note that despite higher SNR, the rms errors for LISA are larger than the ones for TianQin. This is mainly due to the fact that the SMBHB accumulates the majority of its SNR within a few days before the merger (see Fig. 3); hence, it behaves as a relatively short-lived signal, in contrast to galactic white dwarf binaries or EMRIs, in the space-borne detectors. Therefore, the amplitude and polarization modulations (see Eq. 9 and Eq. 16) of the GW signal, in which the information of the source is embedded Tinto and de Araujo (2016), induced by the time varying antenna pattern is more significant for TianQin, since the rotation period of its constellation (3.65 day) is two orders of magnitude shorter than the one for LISA (one year). A comprehensive comparison of the two detectors, in terms of the abilities of detection and parameter estimation for SMBHB mergers, will be analyzed in a followup paper as this is beyond the scope of the current work.

For $z=1.0$ (corresponding to ${D_{L}}\approx 6.8~{}{\rm{Gpc}}$) and $z=2.0$ (corresponding to ${D_{L}}\approx 16~{}{\rm{Gpc}}$), $\Delta\ln M_{c}$ and $\Delta\ln\eta$ can be measured to be $\sim 0.04\%-1.5\%$, $\sim 7.8\%-15\%$, and $\sim 0.1\%-3.5\%$, $\sim 15\%-32\%$, respectively. The sky-position resolution $\Delta\Omega$ can be measured to be $\sim 6.8-94~{}{\deg^{2}}$ and $\sim 35-941~{}{\deg^{2}}$. The error for the estimated time of coalescence $\Delta t_{c}$ is $\sim 15$ sec –$1.0$ hour and $\sim 39$ sec –$4.2$ hour. $\Delta\ln D_{L}$ can be determined to be $\sim 2.3\%-8\%$ and $\sim 5.2\%-24\%$.

The estimation accuracy of the other parameters are also given in Table 1. Moreover, the histograms of the rms errors for $\{\Delta\ln M_{c}$, $\Delta\ln\eta$, $\Delta\phi_{c}$, $\Delta t_{c}$, $\Delta\ln D_{L}$, and $\Delta\Omega\}$ and their fitted distributions for the Monte Carlo simulation of $10^{3}$ trials of ${10^{6}+10^{6}}{M_{\odot}}$ SMBHBs located at $z=0.5$ are shown in Fig. 4 (red lines).

V.2 Effects of spin and eccentricity on parameter estimation

Either spin or eccentricity correction to the phase has been ignored for simplification in many previous works Cutler and Flanagan (1994); Cutler (1998); Poisson and Will (1995); Berti et al. (2005); Arun et al. (2005); Arun (2006). Here, we will show how these effects alter the parameter estimation accuracy for the case of TianQin. Following Lang and Hughes (2008), we adopt a ‘special’ correlation matrix in which diagonal elements are rms errors of parameters while off-diagonal elements are correlation coefficients (see Eq. 29) between each pair of parameters for the 3PN phase (Table 2), the 3PN phase including either eccentricity (Table 3), or spin (Table 4). Different from Lang and Hughes (2008) where the values of matrix elements are given for a specific trial, we present the medians of the rms errors in diagonal elements and the means of the absolute values of correlation coefficients in off-diagonal elements for $10^{3}$ trials of $10^{6}+10^{6}{M_{\odot}}$ SMBHBs. In the simulations, the correlation coefficients change signs in different trials; therefore, the absolute values allow us to take particular note of the magnitude of correlations. We set $z=0.5$ for all sources.

In Table 2, more than half of the correlation coefficients have values $\sim 0.2-0.5$. Three elements have values $\sim 0.7-0.8$, two of which are between $\Delta\ln D_{L}$ with $\Delta\theta$ and $\Delta\phi$. The correlation between mass-related parameters $\Delta\ln M_{c}$ and $\Delta\ln\eta$ turns out to be strongest with coefficient value $>0.9$.

When eccentricity of the binary orbit is considered, the correlation between $\Delta\ln M_{c}$ and $\Delta\ln\eta$ is enhanced somewhat and the rms errors of both $\Delta\ln M_{c}$ and $\Delta\ln\eta$ are approximately doubled (see Table 3). In contrast, the correlation between $\Delta\theta$, $\Delta\phi$, $\Delta D_{L}$ and the other parameters (excluding $\Delta{e_{0}}$) and their rms errors remain almost the same. Hence, the 3D localization error of the source ($\Delta\theta$, $\Delta\phi$ and $\Delta\ln D_{L}$) are not sensitive to the inclusion of eccentricity.

The results for 3PN phase including spin-orbit parameter $\beta$ and spin-spin parameter $\sigma$ are given in Table 4. Among 36 correlation coefficients, 21 have values $\sim 0.2-0.5$, while $\sim 1/6$ have values $0.7-0.9$. Unlike $e_{0}$, $\beta$ and $\sigma$ show strong correlation ($>0.8$) with some of the other parameters, such as $\ln{M_{c}}$, $\ln\eta$, $t_{c}$ and $\phi_{c}$, which worsens the rms error $\Delta\ln{M_{c}}$, $\Delta\ln\eta$, $\Delta{t_{c}}$ and $\Delta{\phi_{c}}$ by a factor of several to several tens ($\Delta\phi_{c}$ increases by a factor of $\sim 70$ considering $c_{\phi_{c},\sigma}=0.986$). As in Table 3, neither the correlation coefficients relevant to $\theta$, $\phi$ and $\ln D_{L}$ nor their rms errors have been changed significantly due to the inclusion of additional parameters. From Fig. 4 one can see that spin (orange curves) has larger impact on the accuracy of parameter estimation than eccentricity (green curves). However, as noted in Vecchio (2004) by using the simple precession model Apostolatos et al. (1994) that the spin induced precession, which is neglected in the current work due to the assumption of the alignment of the spins to the orbital angular momentum, can significantly reduce the rms errors in the parameters, which is a consequence of the new signatures in the GW waveform introduced by precession. The impact of precession on parameter estimation accuracy for TianQin will be the subject of our future investigations.

Comparing the ${10^{6}+10^{6}}{M_{\odot}}$ case ($z=0.5$) in Table 1 with Table 2–4, one can see that $\Delta\ln M_{c}$, $\Delta\ln\eta$, $\Delta t_{c}$, $\Delta\phi_{c}$ in the 3PN phase with eccentricity and spin (3PN+E+S) have been worsened than the ones in 3PN, 3PN+E, and 3PN+S cases. The correlations between mass-related parameters $(M_{c},\eta)$ with $(\phi_{c},t_{c})$ have been discussed in Cornish and Porter (2006); Lang and Hughes (2008); Rodriguez et al. (2013). As found in Lang and Hughes (2008), these parameters will be underestimated if neglecting the correlations among them. However, $\Delta\theta$, $\Delta\phi$, $\Delta D_{L}$ have not been changed significantly. These features can be seen more evidently from Fig. 4.

V.3 Effects of PN order on parameter estimation

The 2PN phase has been adopted as a benchmark in many previous investigations on parameter estimation accuracy of coalescing binaries for either LIGO-Virgo type detectors Poisson and Will (1995); Królak et al. (1995) or LISA Berti et al. (2005); Lang and Hughes (2006). In this part, we will show the difference of the parameter estimation accuracy forecasted by the 2PN and 3PN phases in the context of TianQin. For a comprehensive study, we include eccentricity and spin parameters in the waveforms. The remaining setups of the simulations for the 2PN phase are the same as in Sec. V.1, except $z=0.5$ for all sources in this subsection.

Figure 5 shows the dependencies of the medians of the rms errors on the total mass of SMBHBs when considering the 2PN (green lines with dots) or 3PN (red lines with triangles) phases, respectively. The components are taken from the mass pairs listed in Table 1. We can see that the medians of $\Delta\ln M_{c}$, $\Delta\ln D_{L}$, $\Delta e_{0}$, and $\Delta\Omega$ based on the 2PN phase overlap with the 3PN phase very well in most of the region that covers the lower mass end; whereas discrepancies between the 2PN and 3PN results, within a fraction of the respective values, are shown in these parameters at the higher mass end ($\sim{10^{7}}{M_{\odot}}$). The 2PN phase overestimates the medians of $\Delta\ln\eta$ and $\Delta\ln\mu$ by a factor of 2 around the higher mass end, while it underestimates the medians of $\Delta\phi_{c}$ and $\Delta\sigma$ by a factor of a few around the lower mass end. The 2PN can either underestimate or overestimate $\Delta\beta$ depending on the value of total mass by, at most, a factor of 2. It turns out that $\Delta t_{c}$ is not sensitive to the PN order at all in the concerned total mass range.

VI Conclusions

This work studied the performance on detection and parameter estimation accuracy of SMBHB inspirals for TianQin, a space-borne GW detector working in the millihertz frequency band. By calculating average SNR, we found that TianQin is capable of detecting SMBHBs assembled in the early Universe. The maximum distance reach in terms of cosmological redshift $z>30$. Using the FIM method, we calculated the estimation accuracies of the parameters contained in the GW strain signal model. We adopted the ‘restricted’ 3PN frequency-domain waveform and considered the contributions from spin (spin-orbit and spin-spin) and first-order eccentricity effects to the phase of the GW signal. Using Monte Carlo simulations of $10^{3}$ binaries sampled uniformly in sky location and orientation, we calculated the rms error distributions of parameters and their correlation coefficients for SMBHBs with component masses in the range of $(10^{5},10^{7}){M_{\odot}}$. Furthermore, we studied the effects of eccentricity, spin and PN order on the parameter estimation accuracy for TianQin.

Two sets of codes in Python pyt and Mathematica mat have been developed independently to implement the FIM method. The results from these are consistent with each other within $1\%$ in terms of relative differences of the rms errors. For the sake of computation time, we ignored the terms, when evaluating the elements in Eq. (33), of the derivatives of the amplitude with respect to the coalescence time $t_{c}$, chirp mass $M_{c}$ and symmetric mass ratio $\eta$ considering that they are negligible compared to the derivatives of the phase. To validate this approximation, we used the complete terms including the derivatives of both amplitude and phase to calculate FIM for a subset of the simulations in Mathematica codes, from which the parameter estimation accuracies are consistent with the ones obtained from Mathematica and Python codes that exclude the minor terms.

Simplified assumptions used in this work can be improved in a more careful treatment. Such as, we only took into account the first-order eccentricity contribution to the phase but not to the amplitude. Especially, the latter will become more important when the eccentricity is noticeably larger and enhance significantly the amplitude around the time of periastron. Besides, higher-order eccentricity contribution Moreno-Garrido et al. (1995); Mikóczi (2010); Hinder et al. (2017) to both phase and amplitude may be introduced in the signal model. Moreover, we assumed that the spin parameters are constant during the inspiraling stage of SMBHBs. It may not be true in the real situation when the spin, by chance, is not parallel to the orbital angular momentum; thus, procession may be induced in orbital dynamics Apostolatos et al. (1994); Vecchio (2004). These factors may change the results presented here considerably.

Finally, all the results presented in this work assume a single Michelson-type interferometer configuration. In fact, the space-borne GW detector formed by a nearly equilateral triangular constellation with two-way optic links along each arm is equivalent to two independent Michelson-type interferometers at the frequencies lower than the transfer frequency ($f_{\ast}\approx 0.28$ Hz for TianQin) Vallisneri (2009). For two (or more) interferometers with independent noise, the total SNR is the root of quadratic sum of individual SNRs while the total FIM is the sum of the individual FIMs. We expect that the results from two interferometers will be improved somewhat over the current ones. However, a thorough investigation of this aspect, especially if involving much more sophisticated data combinations from time delay interferometry (TDI) Tinto and Dhurandhar (2005) (which is used to subtract the dominating laser phase noise), is out of the scope of the current paper. All the considerations mentioned above will be subject to our further investigations for TianQin.

Acknowledgements.

Y.W. is supported by the National Natural Science Foundation of China under Grants No. 91636111, No. 11690021, and No. 11503007. Y.M.H. acknowledges the support from the National Natural Science Foundation of China under Grant No. 11703098. This work is partly supported by “the Fundamental Research Funds for the Central Universities” under Grant No. 2019kfyRCPY106. The authors thank the anonymous referee for helpful comments and suggestions.

Appendix A Frequency-domain waveform in SPA

The SPA gives the leading asymptotic behavior of the generalized Fourier integrals in the following form Bender and Orszag (1978)

[TABLE]

where $f(t)$, $g(t)$, $a$, $b$, and $\lambda$ are all real. A point $c\in(a,b)$ is called a stationary point of $g(t)$ if $g^{\prime}(c)=0$. Suppose $f(c)\neq 0$, $g^{\prime}(t)\neq 0$ everywhere else for $t\in(a,b)$ and $g(t)$ is smooth enough to be expanded as a Taylor series. The leading contribution to $I(\lambda)$ comes from a small interval of width $\varepsilon$ surrounding the stationary point $c$ of $g(t)$, such that

[TABLE]

for $\lambda\to+\infty$. To obtain the leading behavior of the integral, we replace $f(t)$ by $f(c)$ and $g(t)$ by $g(c)+{{{g^{(p)}}(c){{(t-c)}^{p}}}/{p!}}$, where ${{g^{(p)}}(c)}\neq 0$ but $g^{\prime}(c)=\cdots={g^{(p-1)}}(c)=0$. Further, we let $s=t-c$ and replace $\varepsilon$ by $\infty$ approximately, then

[TABLE]

To evaluate the integral, we rotate the contour of integration from the real $-s$ axis by an angle $\pi/2p$ if ${g^{(p)}}(c)>0$ and make the substitution

[TABLE]

where $u$ is real. Or rotate the contour by an angle $-\pi/2p$ if $g^{(p)}(c)<0$ and make the substitution

[TABLE]

Thus,

[TABLE]

for $\lambda\to+\infty$. Here the $\Gamma(\cdot)$ is Gamma function. For our case, $p=2$ and $\Gamma(1/2)=\sqrt{\pi}$, so to the leading order

[TABLE]

where we use the term $+\pi/4$ if $g^{\prime\prime}(c)>0$, or $-\pi/4$ if $g^{\prime\prime}(c)<0$.

Given a GW signal $h(t)=A(t)\cos\Phi(t)$, its Fourier transform is

[TABLE]

where

[TABLE]

Using integration by parts for $I_{2}$, we obtain

[TABLE]

where the first term on the right hand side vanishes since $h(t)=0$ when $t\to\pm\infty$. For the second term, we introduce Riemann-Lebesgue lemma: $\int_{a}^{b}f(t)e^{\mathrm{i}\lambda t}dt\to 0$ $(\lambda\to+\infty)$ provided that $\int_{a}^{b}{\left|f(t)\right|dt}$ exists 111This condition is always satisfied since the signal we considered is nonzero only for a finite period of time., such that the second term will also vanish if the following inequality is satisfied

[TABLE]

where $M_{b}$ is the maximal finite boundary value. Note that the constraints $d\ln A/dt\ll d\Phi/dt$ and ${d^{2}}\Phi/{d{t^{2}}}\ll(d\Phi/dt)^{2}$ Creighton and Anderson (2011) are sufficient but not necessary to make the second term vanished. In Appendix B, we will demonstrate that Eq. 44 is satisfied for TianQin. Thus the GW signal in frequency domain is

[TABLE]

Let $\lambda={10^{8}}\times 2\pi f$ in order to use SPA, so the condition $\lambda\to+\infty$ (actually sufficiently large) can always hold in the millihertz frequency band of a space-borne GW detector, such as LISA and TianQin. Suppose $g(t)={10^{-8}}[t-{\Phi(t)}/(2\pi f)]$, then $g^{\prime}(t)={10^{-8}}[1-{\Phi^{\prime}(t)}/(2\pi f)]$, thus the stationary point $t_{\rm{sp}}$ of $g(t)$ is the time at which $d{\Phi(t)}/{dt}=2\pi f$. Furthermore, for the second derivative $g^{\prime\prime}(t)=-{10^{-8}}{\Phi^{\prime\prime}(t)}/(2\pi f)<0$, we use the factor $e^{-i\pi/4}$.

Using Eq. 40, we can obtain $\tilde{h}(f)$ for $f>0$:

[TABLE]

Appendix B Validation of SPA for TianQin

The validity of SPA depends on the satisfaction of Eq. (44) in which the three important terms are the derivatives of amplitude and phase in the time domain with respect to $t$. We set $t_{c}=\phi_{c}=\psi=\iota=0$, $\beta=\sigma=0$, and $e_{0}=0$ in the waveforms. The expressions of the amplitude and the phase of the waveform in the time domain are given by

[TABLE]

with Creighton and Anderson (2011)

[TABLE]

Taking the typical ${10^{6}}+{10^{6}}{M_{\odot}}$ binary system with cosmological redshift $z=0.5$ as an example, we evaluated the left-hand side of Eq. (44) for TianQin and found that its value is well bounded for systems located in $\theta\in(0,\pi)$ and $\phi\in(0,2\pi)$ at the concerned time interval $t\in({t_{\rm{in}}},{t_{\rm{fin}}})$. Here ${t_{\rm{in}}}$ and ${t_{\rm{fin}}}$ are the time that $f={f_{\rm in}}$ and ${f_{\rm fin}}$.