Testing the Kerr nature of supermassive and intermediate-mass black hole binaries using spin-induced multipole moment measurements

The waveform model we employ is the same as that of reference [46] that describes the inspiral dynamics of a non-precessing compact binary system which is based on the post-Newtonian (PN) technique [76–87], and can be schematically represented as,

$$\begin{equation}\tilde {h}\left(f\right)=\frac{{M}^{2}}{{D}_{L}}\sqrt{\frac{5\enspace \pi \enspace \eta }{48}}\sum _{n=0}^{4}{V}_{2}^{n-7/2}\enspace {C}_{2}^{\left(n\right)}\enspace {\mathrm{e}}^{\text{i}\enspace \left({2\enspace {\Psi}}_{\text{SPA}}\left(f/2\right)-\pi /4\right)},\end{equation} \tag{ 2.1 }$$

where Ψ_SPA(f) is the 4PN ⁵ accurate point particle phase which explicitly contain spin-induced quadrupole moment terms at 2PN, 3PN, and 3.5PN and leading order spin-induced octupole moment term at 3.5PN. The leading (second) harmonic and its corrections are incorporated up to 2PN (n = 4) order in the amplitude and appear through the coefficients ${C}_{2}^{\left(n\right)}$ [83]. In equation (2.1), M, η and D_L denote the total mass, the symmetric mass ratio of the binary system and the luminosity distance to the source. The pre-factor V₂ is a function of the total mass of the binary system and the gravitational wave frequency (see section 2 of [46] and supplemental material of [1] for more details about the waveform). Notice that the inclusion of precession effects in the waveform may lead to tighter constraints because of the additional features in the precessing waveforms. However, due to the unavailability of analytical precessing waveform models, we postpone this for future work.

We use Fisher information matrix analysis to obtain the measurement errors on spin-induced quadrupole moment parameters of supermassive and intermediate-mass black holes. Fisher information matrix approach [46, 88, 89] is a semi-analytical parameter estimation technique which can be used to compute the 1 − σ error bars on parameters characterizing the gravitational wave signal, given a waveform model and the detector sensitivity. The elements of the matrix are defined as follows,

$$\begin{equation}{{\Gamma}}_{ij}=2{\int }_{{f}_{\text{lower}}}^{{f}_{\text{upper}}}\enspace \left({\partial }_{i}\tilde {h}\left(f\right)\enspace {\partial }_{j}{\tilde {h}}^{{\ast}}\left(f\right)+{\partial }_{i}{\tilde {h}}^{{\ast}}\left(f\right)\enspace {\partial }_{j}\tilde {h}\left(f\right)\right)\frac{\mathrm{d}f}{{S}_{n}\left(f\right)},\end{equation} \tag{ 2.2 }$$

where we denote the frequency domain gravitational waveform as $\tilde {h}\left(f\right)$ and its partial derivative with respect to the i^th parameter of the binary system as ${\partial }_{i}\tilde {h}\left(f\right)$ (here, ${\partial }_{i}{\tilde {h}}^{{\ast}}\left(f\right)$ is the conjugate of ${\partial }_{i}\tilde {h}\left(f\right)$). In our case the set of parameters in the signal manifold which characterise the compact binary consists of,

$$\begin{equation}\rightarrow {\theta }=\left\{\mathcal{M},\enspace \delta ,\enspace {\chi }_{1},\enspace {\chi }_{2},\enspace {\mathcal{\kappa }}_{s},\enspace {\mathcal{\lambda }}_{s},\enspace {t}_{c},\enspace {\phi }_{c}\right\},\end{equation} \tag{ 2.3 }$$

where t_c , ϕ_c are the time and phase at coalescence and $\mathcal{M}$ and δ are the chirp mass and the asymmetric mass ratio of the system and χ₁, χ₂ are the magnitudes of dimensionless spin parameters. The symmetric combination of the spin-induced quadrupole (κ₁ and κ₂) and octupole moment parameters (λ₁ and λ₂) of the binary system are denoted by κ_s and λ_s , respectively. The 1 − σ error bars on each parameter (equation (2.3)) are computed in the high signal-to-noise ratio (SNR) limit [88, 90] as ${\Delta} \overrightarrow {\theta }=\sqrt{{{\Gamma}}_{ii}^{-1}}$, under the assumption that the detector noise is stationary and Gaussian. Notice that while estimating the errors on κ_s and λ_s , we set the anti-symmetric combination to be zero, i.e., κ_a = λ_a = 0.

The bounds we obtained from the Fisher matrix analysis may be seen as a typical order of magnitude estimates and a detailed study based on Bayesian analysis is required to make a more precise quantification of the same. As the Fisher matrix estimates are expected to agree with those from a numerical sampling of the likelihood in the large SNR limit, we have considered systems that have SNRs of the order of hundreds [88, 91–93]. Further, we make sure that the Fisher matrices used in the analysis are not ill-conditioned and discard those which show numerical issues during inversion.

Gravitational wave detector noise is characterised by the noise power spectral density, S_n (f) in equation (2.2). The noise spectral density used for LISA is given by [75],

$$\begin{equation}{S}_{n}\left(f\right)=\frac{20}{3\enspace {L}^{2}}\left(4\enspace {S}_{n}^{\text{acc}}\left(f\right)+2\enspace {S}_{n}^{\text{loc}}+{S}_{n}^{sn}+{S}_{n}^{\text{omn}}\right)\enspace \left[1+{\left(\frac{2\enspace L\enspace f}{0.41c}\right)}^{2}\right]+{S}_{n}^{\text{gal}},\end{equation} \tag{ 2.4 }$$

where,

$$\begin{align*}\hfill {S}_{n}^{\text{acc}}& =\left\{9{\times}1{0}^{-30}+3.24{\times}1{0}^{-28}\left[{\left(\frac{3{\times}1{0}^{-5}\enspace \mathrm{H}\mathrm{z}}{f}\right)}^{10} +{\left(\frac{1{0}^{-4}\enspace \mathrm{H}\mathrm{z}}{f}\right)}^{2}\right]\right\} \frac{1}{{\left(2\enspace \pi \enspace f\right)}^{4}}\enspace {\mathrm{m}}^{2}\enspace {\mathrm{H}\mathrm{z}}^{-1},\hfill \\ \hfill {S}_{n}^{\text{loc}}& =2.89{\times}1{0}^{-24}\enspace {\mathrm{m}}^{2}\enspace {\mathrm{H}\mathrm{z}}^{-1},\hfill \\ \hfill {S}_{n}^{sn}& =7.92{\times}1{0}^{-23}\enspace {\mathrm{m}}^{2}\enspace {\mathrm{H}\mathrm{z}}^{-1},\hfill \\ \hfill {S}_{n}^{\text{omn}}& =4.00{\times}1{0}^{-24}\enspace {\mathrm{m}}^{2}\enspace {\mathrm{H}\mathrm{z}}^{-1},\hfill \\ \hfill {S}_{n}^{\text{gal}}& =1.633{\times}1{0}^{-44}\enspace {\left(\frac{f}{1\enspace \mathrm{H}\mathrm{z}}\right)}^{-7/3}\mathrm{exp}\left(-{\left(\frac{f}{1.426\enspace \mathrm{m}\mathrm{H}\mathrm{z}}\right)}^{1.183}\right)\hfill \\ \hfill & \enspace \quad {\times}\left(1+\mathrm{tanh}\left(-\frac{f-2.412\enspace \mathrm{m}\mathrm{H}\mathrm{z}}{4.835\enspace \mathrm{m}\mathrm{H}\mathrm{z}}\right)\right)\enspace {\mathrm{H}\mathrm{z}}^{-\mathrm{1}}.\hfill \end{align*}$$

Here ${S}_{n}^{\text{acc}}$, ${S}_{n}^{\text{loc}}$, ${S}_{n}^{sn}$, ${S}_{n}^{\text{omn}}$ and ${S}_{n}^{\text{gal}}$ are due to low-frequency acceleration, local interferometer noise, shot noise, other measurement noise, and galactic confusion noise, respectively. The detector arm length L is fixed to be 2.5 × 10⁹meters and c is the speed of light in meters per second. Noise spectral densities used for DECIGO [73] and DECIGO-B [65] are,

$$\begin{align}\hfill {S}_{n}\left(f\right)=& 7.05{\times}1{0}^{-48}\left(1+{\left(\frac{f}{7.36\enspace \mathrm{H}\mathrm{z}}\right)}^{2}\right)+4.8{\times}1{0}^{-51}\enspace {\left(\frac{f}{1\enspace \mathrm{H}\mathrm{z}}\right)}^{-4}\hfill \\ \hfill & {\times}{\left(1+{\left(\frac{f}{7.36\enspace \mathrm{H}\mathrm{z}}\right)}^{2}\right)}^{-1}+5.33{\times}1{0}^{-52}\enspace {\left(\frac{f}{1\enspace \mathrm{H}\mathrm{z}}\right)}^{-4}\enspace {\mathrm{H}\mathrm{z}}^{-\mathrm{1}},\hfill \end{align} \tag{ 2.5 }$$

and

$$\begin{equation}{S}_{n}\left(f\right)=3.03{\times}1{0}^{-46}\left(1+1.584{\times}1{0}^{-2}{\left(\frac{f}{1\enspace \mathrm{H}\mathrm{z}}\right)}^{-4}+1.584{\times}1{0}^{-3}{\left(\frac{f}{1\enspace \mathrm{H}\mathrm{z}}\right)}^{2}\right)\enspace {\mathrm{H}\mathrm{z}}^{-\mathrm{1}} ,\end{equation} \tag{ 2.6 }$$

respectively.

The lower and upper cut-off frequencies (see equation (2.2)) for the analysis are fixed using the following relations,

$$\begin{align}\hfill {f}_{\text{upper}}& =\mathrm{min}\left({f}_{\mathrm{max}},\enspace {f}_{\text{ISCO}}\right)\hfill \\ \hfill {f}_{\text{lower}}& =\mathrm{max}\left({f}_{\mathrm{min}},\enspace {f}_{\mathrm{4}\enspace \mathrm{y}\mathrm{r}}\right).\hfill \end{align} \tag{ 2.7 }$$

For LISA, we fix f_max as 0.1 Hz and f_min to be 10⁻⁴ Hz and for DECIGO f_max is taken to be 10 Hz throughout the study. To examine the effect of the lower cut-off frequency on different DECIGO configurations, we consider two different scenarios, basic DECIGO or DECIGO-B [65, 69] with a conservative low-frequency cut-off of 10⁻¹ Hz (f₁) and DECIGO [73] at its designed sensitivity (which we refer to as DECIGO) with a low-frequency cut-off of 10⁻² Hz (f₂). As expected, due to the improved sensitivity, bounds obtained from DECIGO are much better than DECIGO-B in general. For the current analysis, the waveform model we use has information about the spin-induced multipole moment parameters only in the inspiral as described in section 2.1. By truncating the analysis at f_upper, given in equation (2.7), we avoid any systematic biases that might arise due to the presence of merger-ringdown phases of the dynamics in the waveform. To obtain the value of f_ISCO, which corresponds to the innermost stable circular orbit (ISCO) frequency of Kerr BHs, we use a fitting formula which is a function of the masses and spins of the binary constituents [1, 46, 94, 95]. Our choice of lower cut-off frequency accounts for the fact that the compact binary system spends four years in each of the detector bands. In order to achieve this, we take ${f}_{\mathrm{4}\enspace \mathrm{y}\mathrm{r}}={f}_{\text{upper}}{\left(1+\frac{{m}_{1}\enspace {m}_{2}}{{\left({m}_{1}+{m}_{2}\right)}^{3}}\enspace 6.6{\times}1{0}^{4}\enspace {\mathrm{T}}^{-\mathrm{1}}\right)}^{-\frac{3}{8}}$, where T is fixed to be 4 yr [96].

Left panel of figure 1 shows the noise PSDs of LISA (black), DECIGO (orange) and DECIGO-B (red) configurations, whereas the right panel shows the variation of signal-to-noise ratios (SNR) with the total mass of the compact binaries in the respective frequency bands. The signal-to-noise ratios are calculated by assuming binary systems optimally oriented at a luminosity distance of 1 Gpc with a mass ratio of 1.1 and dimensionless component spins (0.6, 0.3). Notice the improvement in the signal-to-noise ratios at the high mass end when we choose lower cut-off frequency a factor less (f_low = 10⁻² Hz, red solid curve) than the conservative value (f_low = 10⁻¹ Hz, red dashed curve) in the case of DECIGO configurations.

Zoom In Zoom Out Reset image size

In figure 2, we show the errors on spin-induced quadrupole (solid curve) and octupole moment (dashed curve) parameters as a function of the total mass of the system, where we fix the mass ratio $\left({m}_{1}/{m}_{2}\right)$ to be 1.1 and the dimensionless component spins (χ₁, χ₂) to be (0.6, 0.3). Black curves in figure 2 correspond to supermassive black holes (which LISA is more sensitive to) and, orange and red curves show the errors corresponding to intermediate-mass black holes (which DECIGO detectors will be more sensitive to). As we can see from the figure, errors on both the parameters decrease as a function of total mass initially (irrespective of the detector configuration assumed) and then increase. This is because of the combined effect of SNR and the inspiral truncation frequency of the analysis. As the total mass increases, the signal strength of a binary system with fixed location and orientation in the sky increases, but the upper cut-off frequency decreases as it is inversely related to the total mass, decreasing the number of cycles in the detector band.

Zoom In Zoom Out Reset image size

The results shown in figure 2 are not enough to completely assess the capabilities of LISA and DECIGO to carry out the test of Kerr nature of the binary black hole system, as they correspond to certain representative binary configurations. We repeat the analysis for a simulated population of binary black holes, that may be detected by LISA and DECIGO, and obtain the fraction of the total population this test would yield good constraints for.

The simulated population in our case assumes that the binary black hole merger rate per redshift bin in the observer frame follows the relation,

$$\begin{equation}\frac{\mathrm{d}\mathrm{R}\left(z\right)}{\mathrm{d}z}=\mathcal{R}\left(z\right)\frac{\mathrm{d}{\mathrm{V}}_{\text{c}}\left(z\right)}{dz},\end{equation} \tag{ 3.1 }$$

where R(z) is the number of binary coalescence per observation time, $\mathcal{R}\left(z\right)$ is the merger rate density in the detector frame and V_c(z) denotes the comoving volume.

For LISA sources, we assume the massive black hole rate evolution follows the models given in Klein et al [97]. These semi-analytical massive black hole-galaxy coevolution models assume two different birth mechanisms for massive black holes and also account for the time delay between massive black hole merger and Galaxy merger [98–105]. Following the R(z) and detector frame total mass distributions given in figure [3] of reference [97], we populate binary black holes keeping the component masses nearly equal for three different formation mechanisms described in [97]. ⁶ Among the three models, Model Q3-d (model-1), Model Q3-nod (model-2) and Model popIII (model-3), we observe that very few sources cross the detection threshold (SNR ⩾ 200) for Model popIII and hence we only show the results obtained from Model Q3-d and Model Q3-nod here.

In order to populate intermediate-mass black holes, which are interesting sources for the DECIGO configurations, we start with the following relation [8],

$$\begin{equation}\mathcal{R}\left(z\right)={\mathcal{R}}_{0}\enspace {\left(1+z\right)}^{\lambda }\enspace \frac{{\mathrm{T}}_{\text{obs}}}{1+z},\end{equation} \tag{ 3.2 }$$

where $\mathcal{R}\left(z\right)$ is the rate of binary mergers that occur over the total observation time T_obs measured in the detector frame. Here, ${\mathcal{R}}_{0}$ is a constant which gives the rate density corresponding to a particular value of redshift and we fix it to be 40 Gpc⁻³ yr⁻¹. The magnitude of ${\mathcal{R}}_{0}$ will not affect our results as this is merely a scaling factor. We distribute sources up to a redshift of 20 assuming two different population models for DECIGO configurations, model-1 and model-2. For model-1, we fix λ = 0 (rate density is assumed to be a constant with respect to the redshift) and the component masses to be uniformly distributed in the range 10² − 10⁴ M_⊙. For model-2, we fix λ = 6.5 [8] and the primary mass (m₁) is drawn from a power-law distribution with index 1.6, in the range 10² − 10⁴ M_⊙ and secondary mass (m₂) is drawn from a uniform distribution. The dimensionless spin parameters χ₁ and χ₂ are distributed uniformly between −1 to 1 for both LISA and DECIGO sources. Notice that, among the total populated sources positioned isotropically in orientation and polarisation sky, we choose only those sources which satisfy detection criteria set by the signal-to-noise ratios 200 and 100 for LISA and DECIGO/DECIGO-B respectively. We further perform PE on the signals which pass the signal-noise-ratio threshold using Fisher matrix analysis and the results are shown in figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. Constraints on the BBH nature from spin-induced quadrupole moment measurements

3.2.2. Constraints on the BBH nature from spin-induced quadrupole and octupole moment parameters