Slowly decaying ringdown of a rapidly spinning black hole II:
Inferring the masses and spins of supermassive black holes with LISA

By virtue of the separability of the Hamiltonian for particle’s motion in the Kerr geometry [37], the trajectory of the particle plunging into a Kerr BH while restricted on the equatorial plane ($\theta=\pi/2$) is governed by the following equations:

	$$\displaystyle r^{2}\frac{\,\mathrm{d}t}{\,\mathrm{d}\tau}=-a(a-L_{\mathrm{p}})+\frac{r^{2}+a^{2}}{\Delta(r)}P(r)\>,$$		(II.4)
	$$\displaystyle r^{2}\frac{\,\mathrm{d}\phi}{\,\mathrm{d}\tau}=-(a-L_{\mathrm{p}})+\frac{a}{\Delta(r)}P(r)\>,$$		(II.5)
	$$\displaystyle r^{2}\frac{\,\mathrm{d}r}{\,\mathrm{d}\tau}=-\sqrt{Q(r)}\>,$$		(II.6)

where

	$$\displaystyle P(r):=r^{2}+a^{2}-L_{\mathrm{p}}a\>,$$		(II.7)
	$$\displaystyle Q(r):=2Mr^{3}-{L_{\mathrm{p}}}^{2}r^{2}+2Mr(L_{\mathrm{p}}-a)^{2}\>,$$		(II.8)

$\tau$ is the proper time of the particle and $L_{\mathrm{p}}$ is the orbital angular momentum of the particle normalized by its rest mass $\mu$. We here have taken the energy of the particle normalized by its rest mass $\mu$, $E_{\mathrm{p}}=1$ ²²2This corresponds to an approximation that the particle is at rest at infinity. We are interested in situations where a compact object, originally within the sphere of influence of the SMBH ($\lesssim$ pc for $M\sim 10^{7}~{}M_{\odot}$) but still far away, eventually plunges into the SMBH. The compact object sources ringdown GWs at the moment when it passes the light ring of the SMBH, which is also the case for the particle plunging into the hole from infinity with $E_{\rm p}=1$. Thus, our results shown below are expected to be qualitatively the same, as the trajectory at the light ring will nearly be insensitive to the initial velocity far away from the SMBH.. and the Carter constant to be zero, which ensures that the particle is confined to the equatorial plane for the entire plunging orbit. In our calculations, we assume that the plunging object satisfies the infalling condition: $-2M(1+\sqrt{1+j})<L_{\mathrm{p}}<2M(1+\sqrt{1-j})$ and ignore the self-force of the object under the condition of $\mu\ll M$. Therefore, our waveforms presented here do not comprise any nonlinear effects, e.g., the back-reaction of the gravitational radiation to the light object’s motion. In the following, we use the normalization of $2M=1$.

Figure 1 shows the plunging orbits of a test particle with $L_{\mathrm{p}}=1$ (blue) and $L_{\mathrm{p}}=-1$ (orange) for the cases of $j=0.8,0.9,0.95,$ and $0.99$ (from left to right). A particle with $L_{\mathrm{p}}=1$ travels counterclockwise, and the one with $L_{\mathrm{p}}=-1$ first travels clockwise, and afterward, it is inevitably dragged counterclockwise along the SMBH’s rotation.

The waveform induced by a small object traveling from infinity to the SMBH, as discussed in Sec. II, can be numerically calculated based on the Sasaki-Nakamura (SN) formalism [38]. The perturbation variable of gravitational field, $X_{lm}$, follows the SN equation:

$$\left(\frac{\,\mathrm{d}^{2}}{\,\mathrm{d}{r_{\ast}}^{2}}-F_{lm}\frac{\,\mathrm{d}}{\,\mathrm{d}r_{\ast}}-U_{lm}\right)X_{lm}(\omega,r_{\ast})=\tilde{T}_{lm}(\omega,r_{\ast})\>,$$

(II.9)

where $(l,m)$ denotes a label of the spheroidal harmonic mode, $r_{\ast}$ is the tortoise coordinate, $F_{lm}$ and $U_{lm}$ are functions of a radial coordinate whose explicit forms are given in Ref. [38], and $\tilde{T}_{lm}$ is the source term associated with the geodesic motion of a plunging object whose expression is also given in Ref. [36]. Adopting the Green function method, for a large $r_{\ast}$, Eq. (II.9) can be solved as

$$\begin{split}X_{lm}(\omega,r_{\ast})&=X^{\mathrm{(out)}}_{lm}(\omega)\>e^{\mathrm{i}\omega r_{\ast}}\\ &=\int\tilde{T}_{lm}({r^{\prime}_{\ast}},\omega)G({r^{\prime}_{\ast}},r_{\ast},\omega)\,\mathrm{d}{r^{\prime}_{\ast}}\>,\end{split}$$

(II.10)

where $G({r^{\prime}_{\ast}},r_{\ast},\omega)$ is the Green’s function constructed by the homogeneous solutions of the SN equation. The observed GW signal is expressed by a superposition of all multipole modes,

	$$\displaystyle\tilde{h}(\omega)=\sum_{l,m}\tilde{h}_{lm}(\omega)\>,$$		(II.11)
	$$\displaystyle\tilde{h}_{lm}(\omega):=-\frac{2}{\omega^{2}}{}_{-2}S_{lm}(a\omega,\pi/2)R_{lm}(\omega)\>,$$		(II.12)
	$$\displaystyle R_{lm}(\omega):=\frac{-4\omega^{2}X^{\mathrm{(out)}}_{lm}}{\lambda_{lm}(\lambda_{lm}+2)-6\mathrm{i}\omega-12a^{2}\omega^{2}}\>,$$		(II.13)

where ${}_{-2}S_{lm}$ is the spin-weighted spheroidal harmonics and $\lambda_{lm}$ is the eigenvalue determined by the regularity of ${}_{-2}S_{lm}$. The time domain waveform $h(t)$ can be obtained by the inverse Laplace transform of $\tilde{h}(\omega)$.

Figure 2 shows the numerically computed waveforms for $j=0.99$ and $L_{\rm p}=\pm 1$, decomposed to multipole modes from $(l,m)=(2,2)$ to $(5,5)$. Waveforms for other spin parameters are shown in Appendix. C. From those waveforms, one can see that the higher angular modes are excited more efficiently for rapidly spinning SMBHs, especially for $L_{\mathrm{p}}=1$. As we show in Sec. IV.3, this results in higher angular modes having larger SNRs than the $(2,2)$ mode for rapid spins.

For each multipole mode, we fit the QNM model to the numerical waveform in the time domain by using the least squares fit. The QNM model $h^{\mathrm{(QNM)}}_{lm}$ is given by

$$h^{\mathrm{(QNM)}}_{lm}(t)=\sum_{n=0}^{n_{\mathrm{max}}}A_{lmn}\mathrm{e}^{-\mathrm{i}\omega_{lmn}(t-t_{\ast})+\mathrm{i}\phi_{lmn}}\theta(t-t_{\ast})\>,$$

(III.1)

where $n_{\mathrm{max}}$ is the maximum overtone number in our model, $t_{\ast}$ is the start time of the ringdown, $\omega_{lmn}$ is the QNM frequency for $(l,m,n)$ mode, $A_{lmn}$ and $\phi_{lmn}$ are real-valued amplitude and phase of the mode at $t=t_{\ast}$, and $\theta(x)$ is the step function. Here, $\{A_{lmn},\phi_{lmn},t_{*}\}$ are the fitting parameters for our numerically calculated waveform. Throughout this study, we use qnm, a python module to calculate the QNM frequencies $\omega_{lmn}$ [39].

For each multipole mode of $(l,m)$, we adopt

$$\{(l,m,n)\}=\{(l,m,0),(l,m,1),(l,-m,0),(l,-m,1)\}.$$

(III.2)

To avoid overfitting due to the inclusion of many QNMs in the model, we consider $n_{\rm max}=1$ in this study. The $(l,-m,n)$ modes are called mirror modes, and the excitations of these modes are known to be suppressed for mergers of BHs with comparable masses  [40, 41]. We nevertheless took them into account, as it is less clear whether they are excited for mass ratios of our interest.

For each multiple mode, we assess the goodness-of-fit by evaluating the mismatch ${\cal M}_{lm}$ defined by

$$\mathcal{M}_{lm}(t_{\ast}):=1-\frac{|(h_{lm}(t)|h^{\mathrm{(QNM)}}_{lm}(t))|}{\sqrt{(h_{lm}(t)|h_{lm}(t))(h^{\mathrm{(QNM)}}_{lm}(t)|h^{\mathrm{(QNM)}}_{lm}(t))}}\>.$$

(III.3)

The inner product $(A(t)|B(t))$ is defined by

$$(A(t)|B(t)):=\int^{\infty}_{t_{\ast}}A^{\ast}(t)B(t)\,\mathrm{d}t,$$

(III.4)

which depends on the assumed start time of ringdown $t_{\ast}$. The superscript ^∗ indicates the complex conjugate. For a given $t_{\ast}$, we obtain the best-fit QNM parameters $\{A_{\rm lmn},\phi_{lmn}\}$ as well as the mismatch $\mathcal{M}$. If the amplitude $A_{lmn}$ is extracted self-consistently, $A_{lmn}(t_{*})$ should follow an exponential decay

$$A_{lmn}(t_{\ast})\propto\mathrm{exp}\left[-\frac{t_{\ast}}{\tau_{lmn}}\right]\>,$$

(III.5)

where $\tau_{lmn}:=\{-\mathrm{Im}(\omega_{lmn})\}^{-1}$ is the damping time for the $(l,m,n)$ mode. Comparing $A_{lmn}(t_{*})$ to this relation enables a cross-check that our QNM fitting is valid at sufficiently late times when the relevant QNMs in the fit dominate the ringdown signal.

We here perform the QNM fit for the same cases as in Sec. II. We also consider the case of $L_{\mathrm{p}}=\pm 1$ for each spin $j$ to check the impact of the choice of $L_{\mathrm{p}}$ on the error estimates.

Figure 3 shows the values of mismatch with respect to the assumed start time of ringdown $t_{*}$. One can see that the mismatch decays earlier for lower multipole modes. It is consistent with the proposals [26, 27] that more overtones are excited for higher multipole modes based on the computation of the excitation factors. In all cases, the higher angular modes have better fits than the $(2,2)$ mode due to their larger quality factors, i.e., larger values of $|\text{Re}(\omega_{lmn})/\text{Im}(\omega_{lmn})|$ that enable more robust extraction of individual QNM from the numerical waveforms.

Figure 4 shows the best-fit amplitudes with respect to $t=t_{\ast}$ for the case of $j=0.99$. The figures for other spin cases are shown in Appendix. C. In all cases, the extracted amplitude for the fundamental modes ($n=0$) follows Eq. (III.5), which implies that the amplitude is extracted in a self-consistent manner. For the higher spin cases, the first overtone of the higher multipole modes also follows Eq. (III.5). These fitting results indicate that a large quality factor helps to stably extract the QNM amplitudes. On the other hand, the stable extractions of the mirror modes are not achieved for all cases with $L_{\mathrm{p}}=1$ considered here. However, for the cases with $L_{\mathrm{p}}=-1$ the fundamental mirror modes for higher multipole modes are barely extracted at late times (se Fig. 4). In addition, as shown in Figs. 9, 10, and 11 in Appendix C, the mirror modes are barely extracted for several cases with negative $L_{\mathrm{p}}$, such as the $(4,4)$ and $(5,5)$ modes in the case of $j=0.8$. These facts indicate that infalling objects with negative $L_{\mathrm{p}}$ probably excite the mirror modes to some extent.

The Fisher matrix $\Gamma$ is given by

$$(\Gamma)_{ij}:=\langle\partial_{i}\tilde{h}(f;\bm{\theta})|\partial_{j}\tilde{h}(f;\bm{\theta})\rangle\>,$$

(IV.1)

where $\{\bm{\theta}\}$ is a parameter set in a GW model, $\tilde{h}(f;\bm{\theta})$ is a GW model in frequency domain, and $\partial_{i}$ is the partial derivative with respect to parameter $\theta_{i}$. The inner product $\langle A(f)|B(f)\rangle$ is defined by

$$\langle A(f)|B(f)\rangle:=4\>\mathrm{Re}\left[\int^{f_{\mathrm{max}}}_{f_{\mathrm{min}}}\frac{A^{\ast}(f)B(f)}{S(f)}\,\mathrm{d}f\right]\>,$$

(IV.2)

where $A(f)$ and $B(f)$ are functions of $f$. $S(f)$ is the full noise curve for LISA [47] averaged over the sky position and polarization (the concrete form shown in Appendix. A). We set $f_{\mathrm{min}}=10^{-4}$ Hz which is a conservative choice for the LISA operation, and we take $f_{\mathrm{max}}=1$ Hz, which is large enough for the convergence of the integration in the frequency domain. The RMS errors of $\theta_{i}$, denoted by $\delta\theta_{i}$, and the correlation coefficient between $\theta_{i}$ and $\theta_{j}$, denoted by $c_{\theta_{i}\theta_{j}}$, are estimated by

	$$\displaystyle\delta\theta_{i}=\sqrt{(\Gamma^{-1})_{ii}}\>,$$		(IV.3)
	$$\displaystyle c_{\theta_{i}\theta_{j}}=\frac{(\Gamma^{-1})_{ij}}{\sqrt{(\Gamma^{-1})_{ii}(\Gamma^{-1})_{jj}}}\>,$$		(IV.4)

where $\Gamma^{-1}$ is the inverse matrix of $\Gamma$.

We do not include the mirror modes here because they are subdominant and not extracted stably for all cases as shown in Sec. III.2. Thus, our QNM model for $(l,m)$ mode in the time domain is given by

$$\begin{split}h^{\mathrm{(QNM)}}_{lm}(t)=\>&\big{(}A_{lm0}\mathrm{e}^{-\mathrm{i}\omega_{lm0}(t-t_{\ast})+\mathrm{i}\phi_{lm0}}\\ &+A_{lm1}\mathrm{e}^{-\mathrm{i}\omega_{lm1}(t-t_{\ast})+\mathrm{i}\phi_{lm1}}\big{)}\theta(t-t_{\ast})\>.\end{split}$$

(IV.5)

The QNM model for $(l,m)$ mode in the frequency domain is derived by the Fourier transformation,

$$\tilde{h}^{\mathrm{(QNM)}}_{lm}(\omega)=\frac{\mathrm{i}}{2\pi}\left(\frac{A_{lm0}\mathrm{e}^{\mathrm{i}\omega t_{\ast}+\mathrm{i}\phi_{lm0}}}{\omega-\omega_{lm0}}+\frac{A_{lm1}\mathrm{e}^{\mathrm{i}\omega t_{\ast}+\mathrm{i}\phi_{lm1}}}{\omega-\omega_{lm1}}\right)\>.$$

(IV.6)

Therefore, our QNM model is given by

$$\tilde{h}^{\mathrm{(QNM)}}(\omega)=\sum_{lm}\tilde{h}^{\mathrm{(QNM)}}_{lm}(\omega)\>,$$

(IV.7)

which has 18 model parameters in total, denoted by $\{\bm{\theta}\}$,

$$\{\bm{\theta}\}=\{M,\>j,\>\underset{lmn}{\cup}\{\mathcal{A}_{lmn},\>\phi_{lmn}\}\}\>,$$

(IV.8)

where $\mathcal{A}_{lmn}$ is defined by

$$\mathcal{A}_{lmn}:=\frac{qM_{z}}{d_{\mathrm{L}}}A_{lmn}\>,$$

(IV.9)

and $M_{z}:=(1+z)M$ is the redshifted BH mass.

The value of $t_{\ast}$ in Eq. (IV.6) is chosen at the earliest time when $\mathcal{M}_{lm}\leq 0.05$ ($L_{\mathrm{p}}=-1$) or $\mathcal{M}_{lm}\leq 0.01$ ($L_{\mathrm{p}}=1$) is satisfied for all multipole modes. The fiducial values for $A_{lmn}$ and $\phi_{lmn}$ are set to the best-fit values at $t=t_{\ast}$ for each mode. However, we have confirmed that the final RMS errors for $M$ and $j$ are nearly insensitive to the choice of the upper bound on ${\cal M}_{lm}$ as long as it is $\lesssim 0.1$, whose details are discussed in Appendix. B. Note that in the full Bayesian data analysis, the misalignment between the QNM model and the real signal can lead to systematic bias and leads to the deviation between the estimated values and true values. We here assume that the imposed upper bounds on ${\cal M}_{lm}$ are small enough not to cause the bias in the data analysis.

As we have seen in the previous sections, significant excitations of higher multipole modes are expected for an extreme mass-ratio merger involving a rapidly spinning SMBH. We thus calculate the expected errors of $M$ and $j$ by taking into account higher multipole modes as well. As discussed in Sec. VI in Ref. [32], considering the angular-averaged case, the Fisher matrix can be simplified as³³3When calculating $\Gamma$, we convert the QNM model parameters from the geometrical units to physical units by putting back the dimensional constants $G$, $c$ and $2M$.

$$\Gamma=\Gamma_{22}+\Gamma_{33}+\Gamma_{44}+\Gamma_{55}\>,$$

(IV.10)

where $\Gamma_{lm}$ is the Fisher matrix for each $(l,m)$ mode. This approximation is justified by the approximated relation for the scalar product of spheroidal harmonics [45],

$$\int S^{\ast}_{lmn}(\theta,\phi)S_{l^{\prime}m^{\prime}n^{\prime}}(\theta,\phi)\,\mathrm{d}\Omega\simeq\delta_{ll^{\prime}}\delta_{mm^{\prime}}\>,$$

(IV.11)

where $\,\mathrm{d}\Omega$ is the infinitesimal solid angle, and $\delta_{ij}$ is the Kronecker’s delta.

In this section we discuss the limitations and future avenues of our work, that may be improved by future studies.

An important assumption in our model is neglecting the self-force of the lighter object when calculating the GW waveform. Given that the self-force affects the signal (as well as the final mass and spin of the SMBH) on the order of $q$, we cannot reliably predict precisions better than $\mathcal{O}(q)$ with our waveform model. Therefore, our results are reliable if $\delta\theta_{i}\gtrsim q$ is satisfied, but for $\delta\theta_{i}\lesssim q$ we may need to consider systematic uncertainties due to the waveform modeling. In the cases we investigated for $q=10^{-3}$, this only matters for the near-extremal SMBHs at $d_{\rm L}\lesssim$ a few Gpc (see Fig. 5), and is even less important for smaller mass ratios. Treatment to include these effects is left for future study.

The case of $M\sim 10^{7}M\odot$ and $q\sim 10^{-3}$ or $10^{-5}$ corresponds to a merger between an SMBH and an intermediate-mass BH (IMBH). The event rate of SMBH–IMBH mergers is uncertain, but recent N-body simulations suggest a range $0.003$-$0.03\>\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ [50, 51] or $2$-$20~{}\mathrm{yr}^{-1}$ within $z<1$ ($d_{\mathrm{L}}\lesssim 7\>\mathrm{Gpc}$) [31]. The mergers might be related to the formations of EMRI/IMRIs, where the secondary BH undergoes many inspirals with very high eccentricities before plunging into the BH. Detailed estimates of the event rate as a function of SMBH/IMBH masses are beyond the scope of this work, but combining such rate predictions with our work would enable a more realistic prospect of constraining the properties of SMBHs with LISA.

There is another way to parametrize the QNM model. It utilizes the quality factor, $Q_{lmn}:=\pi f_{lmn}\tau_{lmn}$, rather than $M$ and $j$, and the model is given by

$$h^{\mathrm{QNM}}_{lm}(t)=\sum_{lmn}A_{lmn}e^{-\frac{\pi f_{lmn}}{Q_{lmn}}t}e^{-2i\pi f_{lmn}t+i\phi_{lmn}}\>.$$

(IV.13)

The model parameters here are $\{\bm{{\theta}}\}=\underset{lmn}{\cup}\{f_{lmn},Q_{lmn},A_{lmn},\phi_{lmn}\}$. Using this model we can infer the mass and spin by estimating the multiple peaks of $Q_{lmn}$ in the posteriors and identifying possible combinations of the mass and spin (see Sec. VI.C in Ref. [32] for details). While we have focused on the error estimates introduced in Sec. IV.1 following a more straightforward Fisher matrix analysis, it might be important as future work to consider the prospects of this approach as a complementary way to constrain the SMBH properties.